本發(fā)明涉及彈性力學(xué)領(lǐng)域，具體涉及一種基于Galerkin條形傳遞函數(shù)的薄板振動特性分析方法。

背景技術(shù)：

在彈性力學(xué)中，只有極少數(shù)問題能夠給出解析解，大多數(shù)問題需要通過數(shù)值方法來求解，Kirchhoff板振動特性的求解亦如此。在眾多數(shù)值方法中，有限元法最常見，也是最實用的一種方法。有限元法最大的優(yōu)點是不受求解區(qū)域、邊界條件以及材料屬性的限制，可以分析具有復(fù)雜幾何形狀的彈性力學(xué)問題。但是，有限元方法計算量大，計算時間長，對計算機性能要求高。特別是對某些特殊問題，如彈性體的主動控制的實時計算問題，高梯度應(yīng)力分布的求解問題，以及高頻動態(tài)響應(yīng)計算問題等，有限元法求解精度不高。

條形傳遞函數(shù)方法是一種求解二維彈性力學(xué)問題的半解析數(shù)值方法。這種方法的思想類似于有限條法，也是將求解區(qū)域分為若干個條形區(qū)域，稱為條形單元，在條形單元內(nèi)利用多項式和連續(xù)函數(shù)近似橫向和縱向位移，從而得到基于條形單元的整體微分方程，最后利用傳遞函數(shù)方法求解微分方程，得到半解析解。該方法的一個顯著優(yōu)點是，其既具有有限元方法的靈活性，可以分析復(fù)雜形狀的幾何區(qū)域，同時又能給出封閉形式的高精度半解析解。然而，傳統(tǒng)的條形傳遞函數(shù)方法基于Hamilton原理，需要先給出待求問題對應(yīng)的能量泛函。而問題是，并非所有問題都可以很容易地給出其相應(yīng)的能量泛函，如基于微分型非局部本構(gòu)模型的薄板彎曲問題，這使得條形傳遞函數(shù)方法的應(yīng)用受到了限制。然而，Galerkin方法不需要先寫出待研究問題的能量泛函，可以直接對微分方程進行近似求解。

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

本發(fā)明要解決的技術(shù)問題：針對現(xiàn)有技術(shù)的上述問題，提供一種求解精度高、計算過程數(shù)據(jù)存儲量少、計算效率高的基于Galerkin條形傳遞函數(shù)的薄板振動特性分析方法。

為了解決上述技術(shù)問題，本發(fā)明采用的技術(shù)方案為：

一種基于Galerkin條形傳遞函數(shù)的薄板振動特性分析方法，其步驟包括：

1)將薄板上給定的矩形區(qū)域用NE+1條結(jié)線劃分為NE個矩形區(qū)域，每一個矩形區(qū)域為一個條形單元，第j個條形單元包含第j條結(jié)線和第j+1條結(jié)線以及4個節(jié)點；

2)選取形函數(shù)矩陣N(y)，分別根據(jù)選取的形函數(shù)矩陣N(y)計算每一個條形單元的剛度陣質(zhì)量陣m_e以及荷載向量f_e；

3)基于每一個條形單元的剛度陣質(zhì)量陣m_e以及荷載向量f_e組裝總體 運動微分方程，并針對所述總體運動微分方程處理邊界條件；

4)計算所述總體運動微分方程的傳遞矩陣F(s)以及邊界矩陣M_b(s)與N_b(s)；

5)根據(jù)傳遞矩陣F(s)以及邊界矩陣M_b(s)與N_b(s)計算固有頻率ω。

優(yōu)選地，所述步驟1)中第j個條形單元的結(jié)線位移函數(shù)向量如式(1)所示；第j個條形單元的橫向位移函數(shù)可以通過插值表示如式(2)所示；

φ(x,t)＝{w_j θ_j w_j+1 θ_j+1}^T (1)

式(1)中，φ(x,t)為第j個條形單元的結(jié)線位移函數(shù)向量，w_j為第j條結(jié)線的位移，θ_j為第j條結(jié)線的轉(zhuǎn)角，w_j+1為第j+1條結(jié)線的位移，θ_j+1為第j+1條結(jié)線的轉(zhuǎn)角；

w(x,y,t)＝N(y)φ(x,t) (2)

式(2)中，w(x,y,t)為第j個條形單元的橫向位移函數(shù)，φ(x,t)為第j個條形單元的結(jié)線位移函數(shù)向量，N(y)為形函數(shù)矩陣。

優(yōu)選地，所述步驟2)中選取的形函數(shù)矩陣N(y)為標(biāo)準(zhǔn)Euler梁單元的形函數(shù)。

優(yōu)選地，所述步驟2)中剛度陣的計算函數(shù)表達式如式(3)所示；

式(3)中，為條形單元剛度矩陣，y為條形單元的寬度方向坐標(biāo)軸，l為條形單元的寬度，D為薄板的彎曲剛度，N為選取的形函數(shù)矩陣N(y)，v為薄板的泊松比。

優(yōu)選地，所述步驟2)中質(zhì)量陣m_e的計算函數(shù)表達式如式(4)所示；

式(4)中，m_e為條形單元的質(zhì)量陣，l為條形單元的寬度，ρ為薄板的密度，h為為薄板的厚度，N為選取的形函數(shù)矩陣，y為條形單元的寬度方向坐標(biāo)軸。

優(yōu)選地，所述步驟2)中荷載向量f_e的計算函數(shù)表達式如式(5)所示；

式(5)中，f_e為條形單元的荷載向量；為容許函數(shù)，容許函數(shù)為選取的形函數(shù)矩陣N的轉(zhuǎn)置矩陣N^T，為條形單元結(jié)線上的等效剪力，M_y為條形單元結(jié)線上的y方向彎矩，l為條形單元的寬度。

優(yōu)選地，所述步驟3)中組裝的總體運動微分方程如式(6)所示；

式(6)中，為總體剛度陣，分別由條形單元的剛度矩陣組裝而成，M為由條形單元質(zhì)量陣m_e組成的總體質(zhì)量陣，F(xiàn)為由條形單元荷載向量f_e組成的總體荷載向量，Φ(x,t)為由條形單元結(jié)線位移函數(shù)向量組裝而成的總體位移向量，總體位移向量Φ(x,t)的函數(shù)表達式如式(7)所示；

Φ(x,t)＝{w₁(x,t) θ₁(x,t) w₂(x,t) θ₂(x,t) … w_NE(x,t) θ_NE(x,t)}^T (7)

式(7)中，Φ(x,t)為總體位移向量，w_NE(x,t)為第NE條結(jié)線的位移，θ_NE(x,t)為第NE條結(jié)線的轉(zhuǎn)角。

優(yōu)選地，所述步驟3)中基于每一個條形單元的剛度陣質(zhì)量陣m_e以及荷載向量f_e組成的總體運動微分方程處理邊界條件后如式(8)所示；

式(8)中，分別為處理邊界條件后得到的總體剛度陣，為處理邊界條件后得到的總體質(zhì)量陣，為處理邊界條件后得到的總體荷載向量，為處理邊界條件后的總體位移向量。

優(yōu)選地，所述步驟4)的詳細步驟包括：

4.1)采用式(9)所示函數(shù)表達式計算所述總體運動微分方程的傳遞矩陣F(s)；

式(9)中，I為單位陣，D(s)和D₂的函數(shù)表達式如式(10)所示；

式(10)中，分別為處理邊界條件后得到的總體剛度陣，為處理邊界條件后得到的總體質(zhì)量陣，s為Laplace變換系數(shù)；

4.2)判斷薄板的邊界類型，如果邊界類型為簡支邊界條件，則根據(jù)式(11)計算總體運動微分方程的邊界矩陣M_b(s)與N_b(s)；否則如果邊界類型為固支邊界條件，則根據(jù)式(12)計算總體運動微分方程的邊界矩陣M_b(s)與N_b(s)；

式(11)中，M_b(s)與N_b(s)分別為總體運動微分方程的邊界矩陣，為N₁階單位陣，N₁為未知位移個數(shù)。

式(12)中，M_b(s)與N_b(s)分別為總體運動微分方程的邊界矩陣，為N₁階單位陣，N₁為未知位移個數(shù)。

優(yōu)選地，所述步驟5)中計算薄板自由振動的固有頻率ω的函數(shù)表達式如式(13)所示；

det(M_b(iω)e^-0.5aF(iω)+N_b(ⁱω)e^0.5aF(iω))＝0 (13)

式(13)中，M_b(s)與N_b(s)分別為總體運動微分方程的邊界矩陣，F(xiàn)(s)為總體運動微分方程的傳遞矩陣，a為條形單元長度，det為行列式符號，i為虛數(shù)單位。

本發(fā)明基于Galerkin條形傳遞函數(shù)的薄板振動特性分析方法具有下述優(yōu)點：

1、本發(fā)明直接從薄板平衡方程出發(fā)，利用虛功原理，建立條形單元的運動微分方程，不需要給出薄板對應(yīng)的能量泛函，適用性更廣，計算更簡便；

2、本發(fā)明相對于有限元方法，Galerkin條形傳遞函數(shù)方法在空間的一個方向給出了解析解，在整個空間上給出了半解析解，從而可以顯著地提高求解精度和求解效率；

3、本發(fā)明方法解的形式統(tǒng)一，易于利用計算機編程，并且可以用于相對較復(fù)雜的幾何區(qū)域以及邊界條件問題的求解。

4、本發(fā)明直接從薄板平衡方程出發(fā)，利用虛功原理，建立條形單元的運動微分方程，不需要給出薄板對應(yīng)的能量泛函，不僅適用于基于Galerkin條形傳遞函數(shù)的薄板振動特性分析方法，而且同樣也可用于其他難以直接寫出能量泛函的二維數(shù)學(xué)物理問題的求解。

附圖說明

圖1為本發(fā)明實施例方法的基本流程示意圖。

圖2為本發(fā)明實施例中給定的矩形區(qū)域劃分原理示意圖。

圖3為本發(fā)明實施例中某一個條形單元的結(jié)構(gòu)示意圖。

具體實施方式

如圖1所示，本實施例基于Galerkin條形傳遞函數(shù)的薄板振動特性分析方法的步驟包括：

1)將薄板上給定的矩形區(qū)域用NE+1條結(jié)線劃分為NE個矩形區(qū)域，每一個矩形區(qū)域為一個條形單元，第j個條形單元包含第j條結(jié)線和第j+1條結(jié)線以及4個節(jié)點；參見圖2和圖3，本實施例中給定的矩形區(qū)域長度為a、總寬度為b，用NE+1條結(jié)線劃分為NE個矩形區(qū)域后，每一個條形單元的寬度為l，Oxy為條形單元局部坐標(biāo)系；

2)選取形函數(shù)矩陣N(y)，分別根據(jù)選取的形函數(shù)矩陣N(y)計算每一個條形單元的剛度陣質(zhì)量陣m_e以及荷載向量f_e；

3)基于每一個條形單元的剛度陣質(zhì)量陣m_e以及荷載向量f_e組裝總體運動微分方程，并針對總體運動微分方程處理邊界條件；

4)計算總體運動微分方程的傳遞矩陣F(s)以及邊界矩陣M_b(s)與N_b(s)；

5)根據(jù)傳遞矩陣F(s)以及邊界矩陣M_b(s)與N_b(s)計算固有頻率ω。

本實施例中，步驟1)中第j個條形單元的結(jié)線位移函數(shù)向量如式(1)所示；第j個條形單元的橫向位移函數(shù)可以通過插值表示如式(2)所示；

φ(x,t)＝{w_j θ_j w_j+1 θ_j+1}^T (1)

式(1)中，φ(x,t)為第j個條形單元的結(jié)線位移函數(shù)向量，w_j為第j條結(jié)線的位移，θ_j為第j條結(jié)線的轉(zhuǎn)角，w_j+1為第j+1條結(jié)線的位移，θ_j+1為第j+1條結(jié)線的轉(zhuǎn)角；

w(x,y,t)＝N(y)φ(x,t) (2)

式(2)中，w(x,y,t)為第j個條形單元的橫向位移函數(shù)，φ(x,t)為第j個條形單元的結(jié)線位移函數(shù)向量，N(y)為形函數(shù)矩陣。

本實施例中，步驟2)中選取的形函數(shù)矩陣N(y)為標(biāo)準(zhǔn)Euler梁單元的形函數(shù)，其表達式具體為N＝[N₁ N₂ N₃ N₄]。

本實施例中，步驟2)中剛度陣的計算函數(shù)表達式如式(3)所示；

式(3)中，為條形單元剛度矩陣，y為條形單元的寬度方向坐標(biāo)軸，l為條形單元的寬度，D為薄板的彎曲剛度，N為選取的形函數(shù)矩陣N(y)，v為薄板的泊松比。

本實施例中，步驟2)中質(zhì)量陣m_e的計算函數(shù)表達式如式(4)所示；

式(4)中，m_e為條形單元的質(zhì)量陣，l為條形單元的寬度，ρ為薄板的密度，h為為薄板的厚度，N為選取的形函數(shù)矩陣，y為條形單元的寬度方向坐標(biāo)軸。

本實施例中，步驟2)中荷載向量f_e的計算函數(shù)表達式如式(5)所示；

式(5)中，f_e為條形單元的荷載向量，為容許函數(shù)，容許函數(shù)為選取的形函數(shù)矩陣N的轉(zhuǎn)置矩陣N^T，為條形單元結(jié)線上的等效剪力，M_y為條形單元結(jié)線上的y方向彎矩，l為條形單元的寬度。

對于圖3所示的每個條形單元，用彎矩和扭矩表示的動力學(xué)方程如式(5-1)所示；

式(5-1)中，ρ為薄板的密度，h為為薄板的厚度，w為條形單元的位移，M_x為薄板的x方向彎矩，M_y為薄板的y方向彎矩，M_xy為薄板的扭矩，如式(5-2)所示；

式(5-2)中，D為彈性矩陣，其余符號參數(shù)與式(5-1)相同，彈性矩陣D的具體形式如式(5-3)所示；

式(5-3)中，D為彈性矩陣，D為薄板的彎曲剛度，v為薄板的泊松比。

條形單元上下邊界分別作用等效剪力與彎矩如式(5-4)所示；

式(5-4)中，為剪力，為條形單元結(jié)線上的等效剪力，M_y為條形單元結(jié)線上的彎矩，M_xy為條形單元結(jié)線上的扭矩，l為條形單元的寬度。

式(5-1)所示動力學(xué)方程在y方向上的等效積分“弱”形式如式(5-5)所示；

式(5-5)中，為權(quán)函數(shù)，M_x為條形單元結(jié)線上的x方向彎矩，M_y為條形單元結(jié)線上的y方向彎矩，M_xy為條形單元結(jié)線上的扭矩，l為條形單元的寬度。

式(5-5)等于式(5-6)所示函數(shù)表達式；

式(5-6)中，各個字符參數(shù)含義與式(5-5)相同。

將式(5-2)代入式(5-6)，且令條形單元的位移w等于形函數(shù)N(y)和橫向位移函數(shù)的乘積φ(x,t)(即w＝N(y)φ(x,t))，權(quán)函數(shù)為選取的形函數(shù)矩陣N的轉(zhuǎn)置矩陣N^T即可得到式(5-7)；

式(5-7)中，D為薄板的彎曲剛度，v為薄板的泊松比，N為選取的形函數(shù)矩陣，ρ為薄板的密度，h為為薄板的厚度，x為條形單元的長度方向坐標(biāo)軸，y為條形單元的寬度方向坐標(biāo)軸，為容許函數(shù)，容許函數(shù)為選取的形函數(shù)矩陣N的轉(zhuǎn)置矩陣N^T，l為條形單元的寬度，φ為處理邊界條件后的總體位移。

式(5-7)可記為式(5-8)所示函數(shù)表達式；

式(5-8)中，為條形單元的剛度矩陣，m_e為條形單元的質(zhì)量陣，f_e為條形單元的荷載向量，φ為處理邊界條件后的總體位移。式(5-8)中剛度陣的計算函數(shù)表達式如式(3)所示，質(zhì)量陣m_e的計算函數(shù)表達式如式(4)所示，荷載向量f_e的計算函數(shù)表達式如式(5)所示。

本實施例中，本實施例中，步驟3)中組裝的總體運動微分方程時，將單元剛度矩陣組裝成總體剛度陣的過程和有限元法相同，總體運動微分方程如式(6)所示；

式(6)中，K⁽⁴⁾、K⁽²⁾、K⁽⁰⁾為總體剛度陣，分別由條形單元的剛度矩陣組裝而成，M為由條形單元質(zhì)量陣m_e組成的總體質(zhì)量陣，F(xiàn)為由條形單元荷載向量f_e組成的總體荷載向量，Φ(x,t)為由條形單元結(jié)線位移函數(shù)向量組裝而成的總體位移向量，總體位移向量Φ(x,t)的函數(shù)表達式如式(7)所示；

Φ(x,t)＝{w₁(x,t) θ₁(x,t) w₂(x,t) θ₂(x,t) … w_NE(x,t) θ_NE(x,t)}^T (7)

式(7)中，Φ(x,t)為總體位移向量，w_NE(x,t)為第NE條結(jié)線的位移，θ_NE(x,t)為第NE條結(jié)線的轉(zhuǎn)角。

本實施例中，步驟3)中基于每一個條形單元的剛度陣質(zhì)量陣m_e以及荷載向量f_e組成的總體運動微分方程處理邊界條件后如式(8)所示；

式(8)中，分別為處理邊界條件后得到的總體剛度陣，為處理邊界條件后得到的總體質(zhì)量陣，為處理邊界條件后得到的總體荷載向量，為處理邊界條件后的總體位移向量。

本實施例中，步驟4)的詳細步驟包括：

4.1)采用式(9)所示函數(shù)表達式計算總體運動微分方程的傳遞矩陣F(s)；

式(9)中，I為單位陣，D(s)和D₂的函數(shù)表達式如式(10)所示；

式(10)中，分別為處理邊界條件后得到的總體剛度陣，為處理邊界條件后得到的總體質(zhì)量陣，s為Laplace變換系數(shù)；

4.2)判斷薄板的邊界類型，如果邊界類型為簡支邊界條件，則根據(jù)式(11)計算總體運動微分方程的邊界矩陣M_b(s)與N_b(s)；否則如果邊界類型為固支邊界條件，則根據(jù)式(12)計算總體運動微分方程的邊界矩陣M_b(s)與N_b(s)；

式(11)中，M_b(s)與N_b(s)分別為總體運動微分方程的邊界矩陣，為N₁階單位陣，N₁為未知位移個數(shù)。

式(12)中，M_b(s)與N_b(s)分別為總體運動微分方程的邊界矩陣，為N₁階單位陣，N₁為未知位移個數(shù)。

本實施例中，步驟5)中計算薄板自由振動的固有頻率ω的函數(shù)表達式如式(13)所示；

det(M_b(iω)e^-0.5aF(iω)+N_b(iω)e^0.5aF(iω))＝0 (13)

式(13)中，M_b(s)與N_b(s)分別為總體運動微分方程的邊界矩陣，F(xiàn)(s)為總體運動微分方程的傳遞矩陣，a為條形單元長度，det為行列式符號，i為虛數(shù)單位。

如式(8)所示總體運動微分方程取時間的Laplace變化可得式(13-1)；

式(13-1)中，分別為處理邊界條件后得到的總體剛度陣，為處理邊界條件后得到的總體質(zhì)量陣，為處理邊界條件后的總體位移向量的Laplace變換，為處理邊界條件后得到的總體荷載向量的Laplace變換，s為Laplace變換系數(shù)。

對于薄板的自由振動，總體荷載向量的Laplace變換因此可得式(13-2)；

式(13-2)中，為處理邊界條件后的總體位移向量的Laplace變換，D(s)和D₂的函數(shù)表達式可參見式(10)。

定義狀態(tài)向量η(x,s)如式(13-3)所示；

則式(13-2)可寫成式(13-4)所示函數(shù)表達式；

式(13-4)中，η(x,s)為定義的狀態(tài)向量，F(xiàn)(s)為狀態(tài)矩陣，如式(9)所示?；谑?13-4)，即可推導(dǎo)得出左右邊界條件如式(13-5)所示。

M_b(s)η(-0.5a)+N_b(s)η(0.5a)＝0 (13-5)

式(13-5)中，M_b(s)與N_b(s)稱為邊界矩陣，a為條形單元的長度。判斷薄板的邊界類型，如果邊界類型為簡支邊界條件，則根據(jù)式(11)計算總體運動微分方程的邊界矩陣M_b(s)與N_b(s)；否則如果邊界類型為固支邊界條件，則根據(jù)式(12)計算總體運動微分方程的邊界矩陣M_b(s)與N_b(s)。

求解式(13-4)和式(13-5)，即可得到如式(13-6)所示的解；

η(x,s)＝e^xF(s)(M_b(s)e^-0.5aF(s)+N_b(s)e^0.5aF(s))^-1(13-6)

式(13-6)中，η(x,s)為定義的狀態(tài)向量，F(xiàn)(s)為狀態(tài)矩陣，a為條形單元的長度，M_b(s)與N_b(s)稱為邊界矩陣，最終可推導(dǎo)出式(13-6)的特征方程如式(13)所示。在式(13)所示特征方程的基礎(chǔ)上，令s＝iω，i為虛數(shù)單位，則ω即為薄板自由振動的固有頻率。

綜上所述，本實施例將Galerkin方法與條形傳遞函數(shù)方法相結(jié)合，提出了用于薄板振動特性問題的求解的Galerkin條形傳遞函數(shù)方法，為求解薄板振動特性問題提供一種精度高、計算速度快的新方法。

以上所述僅是本發(fā)明的優(yōu)選實施方式，本發(fā)明的保護范圍并不僅局限于上述實施例，凡屬于本發(fā)明思路下的技術(shù)方案均屬于本發(fā)明的保護范圍。應(yīng)當(dāng)指出，對于本技術(shù)領(lǐng)域的普通技術(shù)人員來說，在不脫離本發(fā)明原理前提下的若干改進和潤飾，這些改進和潤飾也應(yīng)視為本發(fā)明的保護范圍。