基于虛擬元方法求解端部受拋物線荷載的懸臂梁問題

馬俊馳, 梁曉坤, 索宇洋

(遼寧師范大學(xué) 數(shù)學(xué)學(xué)院, 遼寧 大連 116029)

懸臂梁是材料力學(xué)研究中的重要內(nèi)容之一, 廣泛應(yīng)用于建筑工程[1]、電子工業(yè)[2]、機械工程[3]等. 關(guān)于懸臂梁的研究一直備受關(guān)注. 2010年, Yao[4]研究了四階邊值問題(即懸臂梁方程), 證明了一類奇異四階兩點邊值問題正解的局部存在性定理. 文獻[5]利用三角剪切變形理論研究了受拋物線荷載作用的各向同性固定梁. 文獻[6]提出了一種新的無網(wǎng)格局部Petrov-Galerkin(MLPG)方法, 利用一種新的測試函數(shù)來解決二維懸臂梁問題. 由于懸臂梁多樣的荷載條件和復(fù)雜的邊界條件, 直接求解其解析解存在一定的困難. 因此, 尋找一種有效的數(shù)值方法來解決此類問題十分必要. 虛擬元方法[7]可以看作是有限元方法對一般多邊形或多面體單元的擴展, 網(wǎng)格處理比較靈活, 虛擬元空間添加了合適的非多項式函數(shù), 在計算剛度矩陣時, 不需要實際計算這些非多項式函數(shù), 而只需使用自由度.虛擬元空間的引入，便于處理復(fù)雜的單元幾何和高階連續(xù)性條件, 有效地解決更多的實際問題.

1 懸臂梁方程

考慮四階邊值問題

(1)

為了求解問題(1), 對于u∈C4[0,1], 假設(shè)φ(x)=f(x,u(x),u′(x),u″(x),u?(x)),

令v(x)=u″(x), 通過變量代換, 問題(1)可轉(zhuǎn)化成兩個二階問題[8]：

(2)

(3)

問題(1)描述了左端固定右端自由的彈性梁的撓度.在力學(xué)中, 此問題稱為懸臂梁方程.要解決問題(1), 只需求解問題(2)和問題(3).

2 虛擬元方法

2.1 連續(xù)性問題與變分形式

為了研究問題(2)和問題(3),首先考慮二維的泊松方程：

(4)

其中,Ω是2上的一個多邊形區(qū)域, 且f(x,y)∈L2(Ω).方程(4)的變分形式為:找到使得

a(u,v)=(f,v), ?v∈V.

(5)

2.2 構(gòu)造虛擬元空間

首先, 對區(qū)域Ω進行剖分.設(shè){Γh}h是將Ω分解為元素K的序列, 設(shè)εh為Γh的邊的集合,h表示Γh中元素的最大直徑.對于Γh內(nèi)的每個多邊形K,xK為K的重心,hK為K的直徑.當(dāng)k≥1時, 定義空間Βk(?K)為

Βk(?K):={v∈C0(?K):v|e∈Pk(e) ?e??K}.

其中,?K表示多邊形K的邊的集合, 其中，Pk(e)表示次數(shù)不超過k的多項式.定義局部虛擬元空間為

Vk(K)={v∈H1(k):v|?K∈Βk(?K),Δv∈Pk-2},

且P-1(K)={0}.定義全局虛擬元空間為

Vh={v∈H1(Ω):?K∈Γhv|?K∈Βk(?K),Δv|k∈Pk-2}.

在虛擬元方法中, 自由度對形函數(shù)的計算至關(guān)重要, 因此, 給出自由度的計算方法.空間Vk(K)的函數(shù)僅僅由以下自由度決定:①vh在頂點處的值; ②對于k>1,vh在每條邊e上的值; ③對于k>1, 在每個單元上:

其中,

為了方便對雙線性形式進行局部離散, 引入下列投影算子.

定義投影算子Π?:Vk(K)→Pk(K),對于所有的vh∈Vk(K), 都滿足下列正交性條件(?pk,?(Π?vh-vh))0,K=0, ?pk∈Pk(K).

定義L2投影算子Π0:Vk(K)→Pk(K), 對于所有的vh∈Vk(K), 都有(pk,Π0vh-vh)0,K=0, ?pk∈Pk(K).

當(dāng)k=1或k=2時, 顯然有Π?=Π0.通過構(gòu)造與問題的雙線性形式有關(guān)的投影算子, 可以方便離散雙線性形式以及計算剛度矩陣.

2.3 虛擬元離散

基于上述空間和投影算子的定義, 下面將雙線性形式進行局部離散.

局部離散雙線性形式為

aK(u,v)=aK(Π?u,Π?v)+aK(u-Π?u,v-Π?v) ?u,v∈VK,k.

選取φ1,…,φNdof為標(biāo)準(zhǔn)基, 則dofi(φj)=δij,i,j=1,2,…,Ndof.

局部剛度矩陣為

其中,

上述的局部離散雙線性形式, 具有以下兩個重要性質(zhì).

引理1[9]對于任意的多邊形K∈Γh, 都有

(1)一致性:對于?vh∈Vh|K, 并且?p∈Pk(K), 有

(2)穩(wěn)定性:存在不依賴h和K的常數(shù)α*和α*, 使得對?vh∈Vh|K, 有

綜上, 得到離散方程:找到uh∈Vh, 使得

ah(uh,vh)=〈fh,vh〉, ?vh∈Vh.

(6)

由一致性與穩(wěn)定性易證式(6)有且僅有唯一解uh.

2.4 誤差分析

基于上述離散方程, 根據(jù)文獻[7], 可得數(shù)值解的收斂性.下面對H1范數(shù)和L2范數(shù)下數(shù)值解的誤差估計進行討論.

定理1對充分小的h, 問題(6)有唯一解u∈Vh,H1范數(shù)下的誤差估計為

‖u-uh‖1≤Chk(‖u‖k+1+|f|k),

其中,C是一個與h無關(guān)的正常數(shù).

(7)

定理2對充分小的h,L2范數(shù)下的誤差估計如下:

‖u-uh‖0≤Chk+1(‖u‖k+1+|f|k),

其中,C是一個與h無關(guān)的正常數(shù).

因為

a(u-uh,ψ-ψI)≤Chk+1‖u‖k+1‖u-uh‖0,

因此,‖u-uh‖0≤Chk+1(‖u‖k+1+|f|k), 結(jié)論得證.

3 數(shù)值實驗

本節(jié)通過對端部受拋物線荷載作用的懸臂梁進行數(shù)值模擬, 驗證了上述理論分析的準(zhǔn)確性.對于上述懸臂梁, 文獻[11]給出了位移場的精確解.應(yīng)用虛擬元方法求解, 首先對網(wǎng)格進行剖分, 網(wǎng)格生成的具體步驟參考文獻[12].這里選擇的剖分區(qū)域為Ω=[0,8]×[-2,2].圖1給出了網(wǎng)格剖分?jǐn)?shù)N分別為50, 100, 200, 400, 800, 1 600的剖分圖.

圖1 不同網(wǎng)格數(shù)對應(yīng)的剖分圖Fig.1 The division diagram corresponding to different mesh numbers

圖2 不同網(wǎng)格數(shù)對應(yīng)的數(shù)值解Fig.2 Numerical solutions corresponding to different mesh numbers

下面討論不同范數(shù)意義下的相對誤差與絕對誤差.

表1給出了不同網(wǎng)格下數(shù)值解的誤差, 其中, Err(L2)和err(L2)分別表示L2范數(shù)下的絕對誤差與相對誤差, Err(H1)和err(H1)分別表示H1范數(shù)下的絕對誤差與相對誤差.

表1 不同網(wǎng)格數(shù)對應(yīng)的相對誤差與絕對誤差Table 1 Relative error and absolute error corresponding to different mesh number

由表格知, 變化網(wǎng)格剖分個數(shù)N=50, 100, 200, 400, 800, 1 600,L2范數(shù)下的相對誤差和H1范數(shù)下的相對誤差呈現(xiàn)明顯收斂的趨勢,L2范數(shù)下的絕對誤差和H1范數(shù)下的絕對誤差更加精確, 從而驗證了虛擬元離散格式的收斂性與有效性.