本發(fā)明屬于三維力學(xué)振動(dòng)分析數(shù)值求解
技術(shù)領(lǐng)域：
，涉及一種基于疊層基函數(shù)的三維力學(xué)模態(tài)仿真模擬方法。
背景技術(shù)：
：電子器件的使用環(huán)境往往十分惡劣，例如在崎嶇上路上運(yùn)輸時(shí)的振動(dòng)、飛機(jī)起飛、坦克高速行進(jìn)、衛(wèi)星和導(dǎo)彈上升階段的重力加速等對電子器件的機(jī)械強(qiáng)度提出了十分嚴(yán)格的要求。機(jī)械性能卻又是電子器件的可靠性和穩(wěn)定性的重要組成部分，這直接影響到器件能否正常工作。因此對電子器件的機(jī)械性能進(jìn)行優(yōu)化設(shè)計(jì)是有必要的，而模態(tài)分析可以獲得電子器件的振動(dòng)特性，是其機(jī)械性能設(shè)計(jì)的重要環(huán)節(jié)，因此模態(tài)分析中高精度的獲得器件的振動(dòng)特性具有極其重要的意義。目前，利用各種力學(xué)計(jì)算方法對電子器件結(jié)構(gòu)模態(tài)進(jìn)行仿真分析時(shí)，都是采用的有限元本征分析方法。有限元分析一般包括，單元?jiǎng)澐帧卧治?、系統(tǒng)綜合、引入條件、求解方程組、后處理等幾個(gè)步驟。在單元分析中針對單元?jiǎng)澐值玫降拿總€(gè)單元，假設(shè)待求函數(shù)的近似式，選擇若干個(gè)結(jié)點(diǎn)，將近似函數(shù)表示成節(jié)點(diǎn)上插值函數(shù)的形式，然后利用單元上的插值函數(shù)，通過物理分析或數(shù)學(xué)分析，建立單元節(jié)點(diǎn)的待求函數(shù)與外界條件之間滿足的方程組。在這一過程中，插值函數(shù)一般是一階、二階、三階甚至更高階數(shù)，目前各種力學(xué)計(jì)算方法都采用的是[有限元法原理簡明教程，110-117頁，作者：廖日東]一文中提到的畫面法來形成各階基函數(shù)，相對于高階基函數(shù)來說，這種方法實(shí)現(xiàn)起來比價(jià)困難，而且不利于后期加速算法加載。當(dāng)幾何機(jī)構(gòu)比較規(guī)則，結(jié)構(gòu)比較簡單采用低階基函數(shù)來進(jìn)行有限元分析，其精度是可以接受的。隨著分析的結(jié)構(gòu)復(fù)雜化，低階基函數(shù)的精度已經(jīng)不能滿足設(shè)計(jì)者的要求，因此需要構(gòu)造高階基函數(shù)來實(shí)現(xiàn)有限元分析。技術(shù)實(shí)現(xiàn)要素：針對上述存在問題或不足，為解決構(gòu)造有限元模態(tài)分析中的高階插值基函數(shù)，從而獲得高精度的數(shù)值模擬結(jié)果；本發(fā)明提供了一種基于疊層基函數(shù)的三維力學(xué)模態(tài)仿真模擬方法。該基于疊層基函數(shù)的三維力學(xué)模態(tài)仿真模擬方法，包括以下步驟：A.將目標(biāo)電子器件結(jié)構(gòu)進(jìn)行建模，引入位移邊界條件或者應(yīng)力邊界條件建立對應(yīng)的幾何結(jié)構(gòu)模型；B.建立電子器件結(jié)構(gòu)的線性彈性力學(xué)基本邊值問題矩陣形式；C.采用四面體網(wǎng)格剖分求解域；D.選擇標(biāo)量疊層基函數(shù)，將位移在所有網(wǎng)格內(nèi)用標(biāo)量疊層基函數(shù)展開，并運(yùn)用標(biāo)準(zhǔn)變分原理得到目標(biāo)電子器件結(jié)構(gòu)的有限元方程；E.引入結(jié)構(gòu)的慣性力，得到結(jié)構(gòu)的自由振動(dòng)有限元廣義本征方程；F.求E步驟所獲得的本征方程，獲得特征值λ和對應(yīng)的特征向量即振幅向量；G.對F步驟獲得的特征值和對應(yīng)特征向量進(jìn)行后處理獲得振動(dòng)模態(tài)頻率和對應(yīng)振動(dòng)振型。進(jìn)一步優(yōu)選，所述步驟D中，基函數(shù)為標(biāo)量疊層基函數(shù)，基函數(shù)根據(jù)所求問題的復(fù)雜度和計(jì)算結(jié)果的精度要求按照如下疊層規(guī)則進(jìn)行確定：Wp={λ1iλ2jλ3kλ4l|i=0,1,...,pj=0,1,...,p-ik=0,1,...,p-i-jl=0,1,...,p-i-j-k}---(1)]]>Dim(Wp)=(p+1)(p+2)(p+3)6---(2)]]>式(1)(2)中，Wp表示所有基函數(shù)的集合，p表示選取基函數(shù)的階數(shù)，Dim(Wp)表示基函數(shù)的個(gè)數(shù)，i,j,k,l表示上標(biāo)。λ1，λ2，λ3，λ4為四面體網(wǎng)格最基本的標(biāo)量基函數(shù)，其表達(dá)式和推導(dǎo)過程是一種公知，這里不再闡述。一般基函數(shù)會(huì)附著在點(diǎn)、邊、面、體上則基函數(shù)集合表示成Wp=Wnp⊕Wep⊕Wfp⊕Wvp---(3)]]>其中n,e,f,v表示點(diǎn)、邊、面、體。附著在點(diǎn)上的基函數(shù)為Wnp={λ1,λ2,λ3,λ4}Dim(Wnp)=4,(p=1)---(4)]]>附著在邊上的基函數(shù)為Wep={e1,2p,e1,3p,e1,4p,e2,3p,e2,4p,e3,4p}Dim(Wep)=6(p-1),(p=2,3)---(5)]]>其中附著在面上的基函數(shù)為Wfp={f1,2,3p,f1,2,4p,f1,3,4p,f2,3,4p}Dim(Wfp)=4(p-2)(p-1)2---(6)]]>其中附著在體上的基函數(shù)為Wvp={v1,2,3,4p}Dim(Wvp)=(p-3)(p-2)(p-1)6---(7)]]>其中進(jìn)一步優(yōu)選，所述步驟G中后處理具體為：運(yùn)用處理特征值λ，處理對應(yīng)特征向量u為位移矢量，f0為振動(dòng)模態(tài)頻率，N為插值基函數(shù)的矩陣形式。與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果：本發(fā)明提出的基于疊層基函數(shù)的三維力學(xué)模態(tài)仿真模擬方法可以快速的構(gòu)建高階基函數(shù)，且獲得高精度的數(shù)值計(jì)算結(jié)果。附圖說明圖1是本發(fā)明基于疊層基函數(shù)的三維力學(xué)模態(tài)仿真模擬方法的流程圖；圖2是實(shí)施例四節(jié)點(diǎn)四面體單元示意圖；圖3是基函數(shù)附著在點(diǎn)上示意圖；圖4是基函數(shù)附著在邊上示意圖；圖5是基函數(shù)附著在面上示意圖；圖6是基函數(shù)附著在體上示意圖。具體實(shí)施方式下面結(jié)合附圖和實(shí)施例來詳細(xì)說明本發(fā)明的技術(shù)方案。參照圖1，一種基于疊層基函數(shù)的三維力學(xué)模態(tài)仿真模擬方法，包括以下步驟：A.將目標(biāo)電子器件結(jié)構(gòu)進(jìn)行建模，引入位移邊界條件或者應(yīng)力邊界條件建立對應(yīng)的幾何結(jié)構(gòu)模型。根據(jù)電子器件的特性，引入位移邊界條件建立對應(yīng)的幾何結(jié)構(gòu)模型仿真整個(gè)結(jié)構(gòu)的振動(dòng)特性。B.建立電子器件結(jié)構(gòu)的線性彈性力學(xué)基本邊值問題矩陣形式。對于空間邊值問題，在結(jié)構(gòu)(彈性體)內(nèi)部我們要考慮靜力學(xué)、幾何學(xué)、物理學(xué)三方面條件，分別建立三套方程；并給定約束或面力的邊界上，建立位移邊界條件或應(yīng)力邊界條件。具體如下：平衡微分方程∂σx∂x+∂τyx∂y+∂τzx∂z+fx=0∂σy∂y+∂τzy∂z+∂τxy∂x+fy=0∂σz∂z+∂τxz∂x+∂τyz∂y+fz=0---(8)]]>幾何方程ϵx=∂u∂x,ϵy=∂v∂y,ϵz=∂w∂zγxy=∂v∂x+∂u∂y,γyz=∂w∂y+∂v∂z,γzx=∂w∂x+∂u∂z---(9)]]>物理方程σx=E1+μϵx+μE(1+μ)(1-2μ)(ϵx+ϵy+ϵz)σy=E1+μϵy+μE(1+μ)(1-2μ)(ϵx+ϵy+ϵz)σz=E1+μϵz+μE(1+μ)(1-2μ)(ϵx+ϵy+ϵz)τxy=E2(1+μ)γxyτyz=E2(1+μ)γyzτzx=E2(1+μ)γzx---(10)]]>位移邊界條件(u)su=u‾,(v)su=v‾,(w)su=w‾---(11)]]>應(yīng)力邊界條件(lσx+mτyx+nτzx)sσ=p‾x(mσy+nτzy+lτxy)sσ=p‾y(nσz+lτxz+mτyz)sσ=p‾z---(12)]]>上述(8)(9)(10)(11)(12)式中，σx,σy,σz,τxy＝τyx,τyz＝τzy,τzx＝τxz表示求解區(qū)域中6個(gè)應(yīng)力分量，εx,εy,εz,γxy,γyz,γzx表示求解區(qū)域中6個(gè)形變應(yīng)力分量，u,v,w表示求解區(qū)域中3個(gè)位移分量。E是求解區(qū)域中結(jié)構(gòu)的楊氏彈性模量，μ是求解區(qū)域中結(jié)構(gòu)的泊松比，Su表示位移邊界面，Sσ表示應(yīng)力邊界面。l＝cos(n′,x),m＝cos(n′,y),n＝cos(n′,z)，表示應(yīng)力邊界面Sσ上的方向余弦，其中n′為應(yīng)力邊界面Sσ的外法線，x,y,z為應(yīng)力邊界面Sσ上三個(gè)方向的坐標(biāo)值。為位移邊界面Su上的位移值，fx,fy,fz為求解區(qū)域內(nèi)結(jié)構(gòu)受到的各個(gè)方向的體力，為應(yīng)力邊界面Sσ上受到的各個(gè)方向的面力，具體推導(dǎo)過程為一種公知過程，這里不再闡述。在笛卡爾坐標(biāo)系oxyz中，我們假設(shè)彈性體由空間區(qū)域Ω定義，任意點(diǎn)(x,y,z)的位移矢量記為u＝[uvw]T(13)體力矢量記為f＝[fxfyfz]T(14)面力矢量記為p‾=p‾xp‾yp‾zT---(15)]]>應(yīng)變張量記為ε＝[εxεxεxγxyγyzγzx]T(16)應(yīng)力張量記為σ＝[σxσxσxτxyτyzτzx]T(17)需要強(qiáng)調(diào)的是，上述體力矢量，面力矢量，應(yīng)變張量，應(yīng)力張量均為空間坐標(biāo)函數(shù)，式中為簡潔起見未明確寫出。令A(yù)=∂∂x00∂∂y0∂∂z0∂∂y0∂∂x∂∂z000∂∂z0∂∂y∂∂x---(18)]]>L=∂∂x000∂∂y000∂∂z∂∂y∂∂x00∂∂z∂∂y∂∂z0∂∂x=AT---(19)]]>D=E(1-μ)(1+μ)(1-2μ)1μ1-μμ1-μ000μ1-μ1μ1-μ000μ1-μμ1-μ1000000(1-2μ)2(1-μ)000000(1-2μ)2(1-μ)000000(1-2μ)2(1-μ)---(20)]]>n=nx00ny0nz0ny0nxnz000nz0nynx---(21)]]>其中nx＝l,ny＝m,nz＝n則平衡微分方程(8)、幾何方程(9)、物理方程(10)位移邊界(11)、應(yīng)力邊界(12)可分別寫成矩陣形式為Aσ+f＝0(22)ε＝Lu(23)σ＝Dε(24)C.采用四面體網(wǎng)格剖分求解域。采用四面體網(wǎng)格剖分求解域是有限元方法中的一種公知過程，因此本步驟不再詳細(xì)描述。剖分后的求解域被人為分割為多個(gè)三維四面體網(wǎng)格，從而將連續(xù)的幾何結(jié)構(gòu)空間轉(zhuǎn)化為離散的網(wǎng)格空間。D.選擇疊層基函數(shù)，將位移在所有網(wǎng)格內(nèi)用標(biāo)量疊層基函數(shù)展開，并運(yùn)用標(biāo)準(zhǔn)變分原理得到電子器件結(jié)構(gòu)的有限元方程。如果B步驟中的位移矢量u是我們的目標(biāo)函數(shù)，那么對于位移矢量u我們用插值基函數(shù)展開如下形式u=Σi=1nNiuiΣi=1nNiviΣi=1nNiwiT=N100N200...Nn000N100N20...0Nn000N100N2...00Nn·u1v1w1u2v2w2...unvnwnT---(27)]]>令Ni=Ni000Ni000Ni---(28)]]>N＝[N1N2…Nn](29)α＝[u1v1w1u2v2w2…unvnwn]T(30)則式(27)可以記為u＝Nα(31)其中Ni為我們要選擇的插值基函數(shù)。如圖2所示四面體單元中i,j,k,l代表四個(gè)頂點(diǎn)的編號(hào)，我們首先得到四個(gè)最基本的基函數(shù)，具體推導(dǎo)過程為一種公知過程，這里不再闡述：λ1=Ni=ai+bix+ciy+diz6V---(32)]]>λ2=Nj=aj+bjx+cjy+djz6V---(33)]]>λ3=Nk=ak+bkx+cky+dkz6V---(34)]]>λ4=Nl=al+blx+cly+dlz6V---(35)]]>式中ai=xjyjzjxkykzkxlylzl---(36)]]>bi=1yjzj1ykzk1ylzl---(37)]]>ci=1xjzj1xkzk1xlzl---(38)]]>di=1xjyj1xkyk1xlyl---(39)]]>V=161xiyizi1xjyjzj1xkykzk1xlylzl---(40)]]>將(36)式、(37)式、(38)式、(39)式中的i,j,k,l輪換，得到aj,ak,al，bj,bk,bl，cj,ck,cl，dj,dk,dl。V為四面體的體積。為了使四面體的體積不為負(fù)值，單元節(jié)點(diǎn)編號(hào)i,j,k,l必須依照一定的順序。在右手坐標(biāo)系中，當(dāng)i,j,k的方向轉(zhuǎn)動(dòng)時(shí)，右手螺旋應(yīng)指向l的方向。對于標(biāo)量疊層高階基函數(shù)的選擇有如下規(guī)則疊層規(guī)則Wp={λ1iλ2jλ3kλ4l|i=0,1,...,pj=0,1,...,p-ik=0,1,...,p-i-jl=0,1,...,p-i-j-k}---(41)]]>Dim(Wp)=(p+1)(p+2)(p+3)6---(42)]]>式(41)(42)中，Wp表示所有基函數(shù)的集合，p表示選取基函數(shù)的階數(shù)，Dim(Wp)表示基函數(shù)的個(gè)數(shù)，i,j,k,l表示上標(biāo)。一般基函數(shù)會(huì)附著在點(diǎn)、邊、面、體上則基函數(shù)集合表示成Wp=Wnp⊕Wep⊕Wfp⊕Wvp---(43)]]>其中n,e,f,v表示點(diǎn)、邊、面、體。如圖3所示附著在點(diǎn)上的基函數(shù)為Wnp={λ1,λ2,λ3,λ4}Dim(Wnp)=4,(p=1)---(44)]]>如圖4所示附著在邊上的基函數(shù)為Wep={e1,2p,e1,3p,e1,4p,e2,3p,e2,4p,e3,4p}Dim(Wep)=6(p-1),(p=2,3)---(45)]]>其中如圖5(p＝3)所示附著在面上的基函數(shù)為Wfp={f1,2,3p,f1,2,4p,f1,3,4p,f2,3,4p}Dim(Wfp)=4(p-2)(p-1)2---(46)]]>其中如圖6(p＝4)所示附著在體上的基函數(shù)為Wvp={v1,2,3,4p}Dim(Wvp)=(p-3)(p-2)(p-1)6---(47)]]>其中定義彈性體的全部形變勢能密度U1為，具體推導(dǎo)為公知過程這里不再闡述U1=12(σxϵx+σyϵy+σzϵz+τxyγxy+τyzγyz+τzxγzx)---(48)]]>用矩陣表示為U1=12ϵTσ---(49)]]>形變勢能可以全部用位移分量來表示。為此，利用物理方程矩陣式(24)和幾何方程矩陣式(23)代入上式得U1=12(Lu)T(DLu)---(50)]]>則得彈性體全部形變勢能為U=∫∫∫Ω12(Lu)T(DLu)dV---(51)]]>若彈性體受體力和面力作用，空間區(qū)域Ω內(nèi)的體力分量為f＝[fxfyfz]T，sσ邊界上的面力分量為則外力在實(shí)際位移上所做的功稱為外力功W=∫∫∫ΩuTfdV+∫∫sσuTp‾dS---(52)]]>由于外力做了功，消耗了外力勢能，因此在發(fā)生實(shí)際位移時(shí)，彈性體的外力勢能V是V=-W=-∫∫∫ΩuTfdV-∫∫sσuTp‾dS---(53)]]>彈性體的形變勢能與外力勢能之和，即為彈性體的總勢能EPEP＝U+V(54)將式(51)、式(31)和式(53)代入式(54)得EP=∫∫∫Ω12(Lu)T(DLu)dV-∫∫∫ΩuTfdV-∫∫sσuTp‾dS=∫∫∫Ω12(LNα)T(DLNα)dV-∫∫∫Ω(Nα)TfdV-∫∫sσ(Nα)Tp‾dS=∫∫∫Ω12αT(NTLTDLN)αdV-∫∫∫ΩαTNTfdV-∫∫sσαTNTp‾dS=12αT[∫∫∫Ω(NTLTDLN)dV]α-αT∫∫∫ΩNTfdV-αT∫∫sσNTp‾dS---(55)]]>根據(jù)最小勢能原理，即變分原理，彈性力學(xué)微分方程定解問題等價(jià)于要求總勢能EP的最小值，則有即(∫∫∫ΩNTLTDLNdV)α-∫∫∫ΩNTfdV-∫∫sσNTp‾dS=0---(56)]]>也即(∫∫∫ΩNTLTDLNdV)α=∫∫∫ΩNTfdV+∫∫sσNTp‾dS---(57)]]>令K＝∫∫∫ΩNTLTDLNdV(58)F=∫∫∫ΩNTfdV+∫∫sσNTp‾dS---(59)]]>則方程也可簡記為Kα＝F(60)其中K為彈性體的剛度矩陣，α為位移向量，F(xiàn)為外載荷量。E.引入結(jié)構(gòu)的慣性力，得到結(jié)構(gòu)的自由振動(dòng)有限元廣義本征方程。當(dāng)研究結(jié)構(gòu)振動(dòng)問題時(shí)，上述E步驟的α位移向量為時(shí)間的函數(shù)，我們重新定義α(t)＝[u1(t)v1(t)w1(t)u2(t)v2(t)w2(t)…un(t)vn(t)wn(t)]T(61)則根據(jù)E步驟討論得到的有限元方程(60)，引入物體的慣性力得到Mα··(t)+Kα(t)=F(t)---(62)]]>其中M＝∫∫∫ΩρNTNdΩ(63)M為質(zhì)量矩陣，ρ為求解區(qū)域物體的密度，為對時(shí)間的二階導(dǎo)數(shù)。當(dāng)物體自由振動(dòng)時(shí)，此時(shí)F(t)＝0方程(62)退化為Mα··(t)+Kα(t)=0---(64)]]>其振動(dòng)形式叫做自由振動(dòng)，該方程有解的形式α(t)=α^ejωt---(65)]]>這是簡諧振動(dòng)的形式，其中ω為常數(shù)，為振幅向量，將其代入式(64)中，有(-ω2Mα^+Kα^)ejωt=0---(66)]]>消去ejωt后，有Kα^=ω2Mα^---(67)]]>該方程有非零解的條件是|K-ω2M|＝0(68)這就是本征方程。F.求E步驟所獲得的本征方程，獲得特征值λ和對應(yīng)的特征向量即振幅向量。求解E步驟得到的廣義本征方程(68)，得到一系列的特征值λ和對應(yīng)的特征向量G.對F步驟獲得的特征值和對應(yīng)特征向量進(jìn)行后處理獲得振動(dòng)模態(tài)頻率和對應(yīng)振動(dòng)振型。對E步驟獲得特征值λ進(jìn)行處理，對應(yīng)的振動(dòng)模態(tài)頻率為f0=λ2π---(69)]]>根據(jù)得到的本征方程(68)的特征向量結(jié)合基函數(shù)，由公式(31)可以得到求解域內(nèi)的位移分布，這就是對應(yīng)振動(dòng)模態(tài)頻率的振動(dòng)振型。當(dāng)前第1頁1 2 3