周期多孔結構的Steklov彈性特征值問題的多尺度漸近分析

譚理琴， 馬 強， 胡 兵

(四川大學數學學院， 成都 610064)

1 引 言

多孔材料在現代工程中應用非常廣泛.多孔材料微觀結構的自由振動行為的預測在數學上可以歸結為一個求解彈性特征值方程的初邊值問題.該問題通常只能進行數值計算. 對于這類高度非均勻材料，傳統的數值方法需要產生精細的計算網格才能準確捕捉結構的微觀特征，計算量大. 為有效減少計算規模，國內外學者們采用多尺度方法來對結構的微觀構型與宏觀性能之間的關系展開研究[1-3]. 1964年，Marchenko和Khruslov[4]首次提出了多尺度漸近展開法，并針對各類數學物理問題建立了多尺度漸近表達式[5]. 之后，Hou和Wu提出了多尺度有限元方法[6]，并展開了實際計算. 崔和曹[7, 8]引入了復合材料物理力學行為的二階雙尺度(Second-Order Two-Scale, SOTS)漸近展開方法[6-8]，可以更準確捕捉結構的局部微觀振蕩行為.

本文采用二階雙尺度漸近展開方法研究周期多孔結構中的Steklov彈性特征值問題. 以往的研究主要討論了周期性橢圓特征值問題的均勻化方法[9，10]，如Vanninathan[11]給出了在周期孔洞區域上Dirichlet和Steklov以及Neumann三種孔洞邊界條件下的特征值問題，其中的結果都未曾考慮到對特征函數與特征值進行二階雙尺度分析. 我們以往的工作[12]表明二階雙尺度漸近展開式對于有效給出結構特征值與特征函數的近似解是必不可少的.本文討論Stkelov彈性特征值問題的二階雙尺度漸近分析方法，給出結構特征值與特征函數二階展開式，并建立有限元計算方法，最后通過數值算例驗證了多尺度模型有效性.

2 多孔結構Steklov特征值問題

考慮圖1所示的周期孔洞區域Ωε(圖1(a))和單胞區域Q*(圖1(b)). 將孔洞填滿, 我們可得到一個單連通有界區域Ω?RN，其中N表示維度(N=2,3)，邊界?Ω滿足Lipschitz條件. 孔洞定義為ωε，Ωε可表示為Ωε=Ω-ωε. Ωε對應的參考單胞為Q*=Q-ω，其中Q=[0,1]N，ω?Q為有界開集. 這樣, Ωε由單胞Q*以周期ε沿坐標軸排開而形成，即

(a) 位移變形圖(a) Deformation graph

(a) Ωε

其中I是指標集，表示為

I={z∈Zn|ε(Q*+z)?Ωε}.

可以看到，ε為結構宏觀尺度與微觀尺度的比值.同樣得到

此外，假定ω的邊界為S且不與外邊界?Q相交，得到

?Q*=?Q∪S,

基于以上設置，考慮如下周期孔洞區域上的Steklov彈性特征值問題：

(1)

ξ=(ξi∈RN,0≤ρ0<ρε(x)≤M，

這里β1,β2和M是正常數. 由結構的周期性, 這些材料系數可表示為

其中cijkl(y),ρ(y)是關于y的周期函數.?Ω由齊次Dirichlet(Γ1)和Neumann(Γ2)構成，且有?Ω=Γ1∪Γ2,Γ1∩Γ2=?.

上述特征值問題對應的變分形式為: 求(λε,uε)，使得

a(uε,v)Ωε=λε(ρε;uε,v)Sε,?v∈V(Ωε)

(2)

其中

V(Ωε)為Sobolev空間

V(Ωε)={v|v∈(H1(Ωε))N,v=0， on Γ1},

并在此空間中定義范數

3 二階雙尺度分析

(4)

(5)

由方程組可以定義一階校正項表達式

(7)

其中

Nα1(y)=(Nα1km(y))1≤k,m≤N∈W(Q*),

W(Q*)={vij∈(H1(Q*))N×N|vij(y)=vij(y+1),

Nα1km(y)滿足單胞問題

(9)

(10)

這里|·|表示區域測度. 進一步，我們得到均勻化Steklov彈性特征值問題為

(11)

至此我們完成了均勻化過程以預測Steklov彈性特征值問題的宏觀行為.

α2=1,2,...,N

(12)

其中

Nα1α2(y)=(Nα1α2km(y))1≤k,m≤N∈W(Q*).

(14)

(15)

(16)

(18)

由于ωε(x)的變分形式可寫成

a(ωε,v)Ωε=ελ0(ρε;u0,v)Sε

(19)

ελ0(ρε;u0,uε)Sε=λε(ρε;uε,ωε)Sε

(20)

該方程被稱作“校正方程”. 將ωε(x)和uε(x),λε的展開式代入校正方程，比較方程兩端ε各冪次系數得到

O(ε0):λ0(ρε;u0,u0)Sε=λ0(ρε;u0,u0)Sε，

O(ε1):λ0(ρε;u0,u1)Sε=λ0(ρε;u0,u1+ω1)Sε+

λ1(ρε;u0,u0)Sε，

O(ε2):λ0(ρε;u0,u2)Sε=λ0(ρε;u0,u2+ω2)Sε+

λ0(ρε;u1,ω1)Sε+λ1(ρε;u0,u1+ω1)Sε+

λ2(ρε;u0,u0)Sε.

(21)

這樣我們有如下結果：

定理3.1令λε為滿足問題的特征值.則存在與ε無關的常數C，使得

(23)

證明 首先,有

當ε足夠小時，由均勻化理論有

uε(x)→u0(x),ωε(x)→u0(x).

因而ε(ρε;uε,ωε)Sε有上界. 再由Cauchy-Schwarz不等式即可得到式(23). 證畢.

定理3.2令λε為滿足問題的特征值.則存在與ε無關的常數C，使得

|ε-1λε-λε,1|≤

ε3/2‖ωε-u0‖L2(ρε;Sε)+

(24)

|ε-1λε-λε,2|≤

|ε-1λε-λε,1|+ε5/2‖ωε-u0‖L2(ρε;Sε)+

(25)

證明 首先引入函數zε∈V(Ωε)為如下問題的唯一解：

‖zε‖L2(ρε;Sε)≤C‖ε-1λεωε-λ0u0‖L2(ρε;Sε).

于是

‖zε+ωε-u0‖L2(ρε;Sε)≤

C‖ωε-u0‖L2(ρε;Sε).

因此，

ε-1λε-λε,1=

對于上面的等式右邊前兩項，我們分別有如下估計

Cε3/2‖ωε-u0‖L2(ρε;Sε)

(26)

對于第三項，我們有

ε3/2‖ωε-u0‖L2(ρε;Sε)]

(27)

聯合不等式(26)和(27)即得估計式(24).同理也可得到式(25), 證畢.

根據估計式(23)～(25)，對特征值的誤差估計轉化為估計ωε與其SOTS展開式之間在L2(ρε;Sε)范數下的誤差.ωε的雙尺度展開是有效的，并且能夠證明在范數‖·‖V(Ω)下嚴格收斂到u0. 假設有如下正則性：

Nα1,Nα1α2∈W(Q*)∩L2(ρ;S),

u0∈(H2(Ωε))N∩L2(ρε;Sε)，

ωε∈(H1(Ωε))N∩L2(ρε;Sε)，

我們可以得到如下估計：

最后得到如下誤差估計

|ε-1λε-λ0|≤Cε,

|ε-1λε-λε,1|≤Cε2,

|ε-1λε-λε,2|≤Cε2.

4 有限元算法

基于上述二階雙尺度漸近展開模型，我們來建立關于Steklov彈性特征值問題的有限元算法.

(ii) 在單胞區域Q*上利用如下變分形式計算一階單胞函數Nα1km(y):

(28)

(iv) 通過下面單胞區域上的變分形式解得二階單胞函數Nα1α2km(y):

(29)

(v) 基于如下變分公式使用子空間迭代法在Ω中求解均勻化Steklov彈性特征值問題：

5 數值算例

對于圖2所示的二維周期孔洞材料，我們考慮無量綱計算，其中非孔洞區域由單一均勻材料構成，材料楊氏模量為E=2.0e6，泊松比為ν=0.3，密度為ρ=1.0，圖3給出了單胞區域Q*. 在外邊界?Ω上令uε=0，我們針對該結構內部孔洞的平面應變Steklov特征值問題進行漸近多尺度計算.

圖2 ε=1/8時宏觀區域ΩεFig.2 The macroscopic domain Ωε with ε=1/8

圖3 單胞區域Q*Fig.3 The composite cell domain Q*

先將單胞區域、均勻化區域與周期復合區域進行三角剖分，表1顯示了ε=1/8和ε=1/16時的網格信息. 在本算例中, 針對單胞問題我們采用分片線性元，均勻化問題采用分片二次元.對于原多孔區域的Steklov問題, 為減小計算開銷我們也采用分片線性元. 可以看到， 當ε減小時，宏觀細網格規模增大，網格尺寸減小，宏觀細網格有限元計算量增加.使用雙尺度有限元方法后，單胞網格和均勻化網格的單元數與結點數均保持不變，計算時間固定，當區域周期ε=1/16時，單胞個數增多，雙尺度有限元法更具計算優勢.

表1 網格信息Tab.1 Mesh information

在單胞區域中計算得到Nα1km(y)，進一步由材料均勻化公式計算得到均勻化密度ρ0=2.67以及均勻化本構矩陣為

設置收斂容許誤差為10-10， 采用子空間迭代法計算ε=1/8與ε=1/16時模型的前6個特征值. 表2與表3列出了特征值的各階近似解以及相對誤差. 可以看到，當特征階數增加時，我們所使用的二階雙尺度有限元方法得到的各階近似解的相對誤差也會變大. 雖然均勻化特征值和一階校正的結果之間相差并不大，但與細網格下的特征值之間差距較大而經過二階校正后所得到的結果

表3 ε=1/16時前6個特征值和相對誤差Tab.3 The first six eigenvalues and relative errors with ε=1/16

明顯好于一階. 當ε變小的時候，所有細網格下的特征值與各階逼近解之間的相對誤差也減小.

數值算例顯示，通過對各特征值近似解進行定量比較，我們可以驗證二階雙尺度模型的收斂性. 由于各階特征函數的不唯一性，這里我們只進行了簡單的定性比較，也可以對求得的各階特征函數近似解進行最小二乘逼近從而得到它們的定量比較.我們相信, 通過額外的計算二階單胞函數，二階近似解會具有更高的精度.

6 結 論

本文發展了周期孔洞結構的Steklov彈性特征值問題的二階雙尺度漸近分析方法，給出了特征值與對應的特征函數的SOTS漸近展開表達式，對特征值進行了誤差估計，構造了相應的有限元計算方法. 數值算例結果驗證了二階雙尺度分析有限元方法解決周期多孔結構模型Steklov彈性特征值問題的有效性. 我們可以得到如下結論：

(i) 均勻化解只能反映出原問題解的宏觀行為，因而有必要通過一階和二階校正解來反映材料的微觀特性，而得到的二階逼近解顯然比一階逼近解更能捕捉到材料的局部振蕩行為，所以將二階雙尺度漸近展開到二階是合理且必要的；

(ii) 數值算例考慮了前6個特征值，當特征值變大時，使用的二階雙尺度有限元算法得到的各階逼近解的相對誤差也都變大，因而研究特征值的漸近性態對尋找適當的高階逼近解具有重要意義；

(iii) 二階雙尺度有限元算法對于單胞問題和均勻化問題的網格劃分相對來說較為粗糙，且計算時間所占內存固定，使得當ε越趨向于0時宏觀區域的網格尺寸變小，計算量與內存占量變大，因此二階雙尺度方法更具計算優勢.