葛志昊, 吳慧麗

(河南大學數(shù)學與統(tǒng)計學院，開封 475004)

0 引言

連續(xù)介質(zhì)力學中的非局部理論是由文獻[1–3]提出的。近場動力學模型和非局部擴散模型都屬于非局部理論的范疇，自Silling 在文獻[3]中首次介紹近場動力學模型之后，其有效性已經(jīng)在許多領域得到證明，如復合材料的斷裂、多晶體的斷裂、納米纖維網(wǎng)絡、裂縫的不穩(wěn)定、圖像處理等[4–6]。非局部模型是一個不涉及位移場空間導數(shù)的積分型方程，因此，非局部模型允許有跳躍間斷點的解存在，并且可以用同一個非局部方程來描述間斷點及非間斷點處的物理過程，而不需要知道間斷點的位置。此外，可以將非局部模型看作分子動力學的變形，它可以從微觀上反映物理材料中粒子之間的相互作用，將連續(xù)力學與分子動力學聯(lián)系起來[7]，克服了經(jīng)典的連續(xù)模型只能從宏觀的角度研究材料的結(jié)構(gòu)性質(zhì)和影響規(guī)律的缺點。對于近場動力學模型Cauchy 問題的適定性研究，可以參考文獻[8–11]。特別是，文獻[10,12]提出了一種非局部向量微積分計算和抽象的非局部平衡法則。此外，也有一些數(shù)值方法求解該問題，包括分片常數(shù)有限元方法、有限差分法、求積法和質(zhì)點法等[13–16]。最近，文獻[17]提出了分片常數(shù)有限元方法來解決非局部擴散問題。為了得到一種高階數(shù)值方法，我們提出了一種高階有限元方法來求解二維體積約束的非局部擴散問題，并采用一種新技巧計算線性元的剛度矩陣。值得一提的是，求解二維體積約束的非局部擴散問題并不是平凡的。

本文第1 部分，我們提出了體積約束的非局部擴散問題的有限元方法，給出了元素的編碼原理和數(shù)值計算節(jié)點的編碼表達式。此外，我們還采用了一種新的方法來計算剛度矩陣。最后，通過數(shù)值算例驗證了理論計算結(jié)果。

1 基于新技巧的有限元方法

令? ?Rd(d= 2,3)表示一個有界連通開區(qū)域，且?=?s ∪?I，其中?s是解區(qū)域，?I是?s周圍寬度為δ的限制區(qū)域。現(xiàn)在，我們列出了一些如下的符號，關于定義的細節(jié)和一些例子，可以參考文獻[12]。

β=β(x,x′):?×? →R 為一個兩點標量函數(shù)，并且滿足：

注意，δ/2 可以用其他小于δ的適當?shù)某?shù)來代替，以保證β(x,x′)的非退化性。

令α=α(x,x′) :?×? →? ≤Rd是一個反對稱兩點向量函數(shù)滿足α(x,x′)+α(x′,x) = 0。給定一個兩點函數(shù)v :?×? →? ?Rd和一個點函數(shù)u:? →R，非局部散度算子D(v):? →R 和其共軛算子D?(u):?×? →? ?Rd，定義如下

本文研究帶有體積約束的非局部擴散問題

其中?=(?δ,1+δ)×(?δ,1+δ),?s=(0,1)×(0,1)。

問題(1)中的核函數(shù)用γ(x,y)來表示，即γ(x,y)=βα·α。若取

則有

將區(qū)域?s剖分成2N1N2(N1=N2=N)個三角形單元，從而得到步長

同時，以相同的步長對帶狀邊界區(qū)域?/?s進行剖分，將區(qū)域?上的三角形剖分單元的集合記為P=∪{K}=P1∪P2，其中K為三角形剖分單元，P1是求解區(qū)域?s上三角形剖分單元的集合，P2是帶狀邊界區(qū)域上?/?s三角形剖分單元的集合。

記?上有限元節(jié)點的指標集為N，?s上的有限元節(jié)點指標集為NI，?/?s上有限元節(jié)點的指標集為Nb，求解區(qū)域?s上的有限元節(jié)點總數(shù)為n= (N1+1)(N2+1)，區(qū)域?上有限元節(jié)點總數(shù)為n1。

設P和T分別是2×n和3×(2N1N2)的矩陣，矩陣P的第k列表示三角形剖分中整體編碼是k的節(jié)點的坐標，矩陣T的第k列表示三角形剖分中編碼是k的三角形剖分單元的三個頂點對應的整體編碼。剖分節(jié)點的整體編碼按照由左至右、由下至上的原則，三角形剖分單元上三個節(jié)點的局部編碼按照逆時針的原則。例如，N1=N2= 2 時，每個三角形剖分單元上剖分節(jié)點的局部編碼如圖1，剖分節(jié)點以及三角形剖分單元的整體編碼如圖2。

圖2 剖分單元及節(jié)點的整體編碼

下面采用線性有限元方法來求解非局部擴散問題(1)，V和V h ?V分別是非局部擴散問題(1)的解空間及有限元子空間。在參考三角形單元上三個頂點對應的節(jié)點局部基函數(shù)分別為

令Am1=(xm1,ym1)′,Am2=(xm2,ym2)′,Am3=(xm3,ym3)′。然后，根據(jù)參考單元與實際單元Km=△Am1Am2Am3之間的仿射變換關系式(2)，我們可以定義實際三角形剖分單元Km上三個節(jié)點對應的局部基函數(shù)ψmi(x,y)=,i=1,2,3，其中

|Jm|=(xm2?xm1)(ym3?ym1)?(xm3?xm1)(ym2?ym1)為Jacobi 矩陣的行列式。因此，整體線性基函數(shù)?i滿足

記?i??i(x)，則非局部擴散問題(1)的變分形式及有限元逼近的一般形式分別是：求u(x)∈V和uh(x)∈V h，使得

其中

uj表示在節(jié)點xj處的數(shù)值解的值。

令v(x)=?i(x),i ∈NI，并且利用式(4)，則有

因此，有

其中C(NI)和C(Nb)分別是NI和Nb索引集的元素的個數(shù)，i,j ∈NI,m ∈Nb。然后，得到線性代數(shù)系統(tǒng)

其中剛度矩陣A 為

在式(9)中，矩陣A1很容易用通常的方法進行計算，矩陣A2可以計算如下

首先，從式(10)取j ∈N 和i ∈NI。然后，可以得到一個n×n1矩陣

其次，從矩陣B中選取包含?j(y)(j ∈Nb)的元素，矩陣B的其余元素保留并表示為B1，而新矩陣B1恰好就是A2(i,j ∈NI)。據(jù)我們所知，這是第一次用B來計算A2。

采用與文獻[18]相似的論證方法，得到如下結(jié)論，此處不再詳細證明。

定理1如果存在常數(shù)s ∈(0,1)，以及γ?和γ?，使得

且u ∈V ∩Ht+1(?)，則對于足夠小的h，存在一個常數(shù)C，使得

2 數(shù)值算例

我們考慮以下二維非局部擴散問題

取?s=(0,1)×(0,1),?=(?δ,1+δ)×(?δ,1+δ)，選擇真解為u(x)=x2+y2，則f(x) =和g(x)可以相應的計算得出。表1給出了L2范數(shù)誤差和H1半范數(shù)誤差，其中N為x(或y)方向等分的份數(shù)。

表1 δ =0.2

從表1 中可以看出，當網(wǎng)格尺寸H變小時，對于分片線性元，L2范數(shù)的誤差收斂階為2，H1范數(shù)的誤差收斂階為1，而文獻[17]中P0元的L2范數(shù)和H1范數(shù)的誤差收斂階都是0.5。因此，可以說本文所設計的方法的收斂階幾乎是最優(yōu)的。從圖3 中我們看到，當網(wǎng)格尺寸h變小時，數(shù)值解更接近真實解。

圖3 δ =0.2 時的精確解和數(shù)值解

3 結(jié)論

本文針對二維體積約束的非局部擴散問題提出了一種基于新技巧的高階有限元方法，該數(shù)值方法的剛度矩陣由一個新的矩陣B提取，該矩陣易于計算。據(jù)我們所知，用矩陣B來計算A2這是第一次，本文給出了元素的編碼原理和數(shù)值計算節(jié)點的編碼表達式。同時，通過數(shù)值算例，驗證了對于體積約束的非局部擴散問題中線性有限元方法的誤差L2范數(shù)的收斂速階達到2，H1半范數(shù)的收斂速階幾乎為1。在今后的工作中，將用本文的數(shù)值方法來處理Pk(k ≥2)元素。