葉舒愉， 馬 強， 胡 兵

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

1 引 言

在物理、力學和工程領域，很多問題在數(shù)學上都被歸結(jié)為求解偏微分方程特征值問題. 對于此類問題，傳統(tǒng)的數(shù)值方法僅考慮對宏觀現(xiàn)象的模擬，難以揭示微觀構(gòu)造的特性，因而對于復雜的材料結(jié)構(gòu)傳統(tǒng)的分析方法不再適用. 于是，國內(nèi)外學者展開了對多尺度方法的研究. 多尺度方法將問題劃分為兩類不同的區(qū)域：第一類是耦合區(qū)域，主要對微-納米構(gòu)件和材料在微-納米尺度下進行研究，主要方法有準連續(xù)介質(zhì)力學方法[1]、原子/連續(xù)耦合方法[2-3]等；第二類是材料區(qū)域，該區(qū)域僅包含原子或連續(xù)介質(zhì)，考慮材料的宏觀性能與細觀結(jié)構(gòu)的關系，代表性方法有均勻化方法，它包括多尺度漸近展開法[4]、收斂估計法[5]、隨機處理法[6]等.

對于多孔材料，由于孔洞區(qū)域結(jié)構(gòu)復雜，本文考慮使用多尺度漸近展開法求解其特征值問題. 1964年，Marchenko和Khruslylov[4]首次提出多尺度漸近展開法. 隨后，Ngnetsent[7]和Allaire[8]提出了雙尺度收斂的概念和雙尺度收斂分析方法，提供了一種獲得均勻化模型的方法. 在此基礎上，Cui和Cao[9]通過對二階展開項的考慮證明了二階雙尺度(Second-Order Two-Scale, SOTS)方法可以有效地預測復合材料的物理力學性能.

此外，使用均勻化方法研究特征值也是一個研究熱點.如Kesavan[10-11]提出了周期情況下橢圓型特征值問題的求解方法，Vanninathan[12]研究了周期性孔洞區(qū)域的特征值問題，Bonnetier等[13]討論了Neumann-Poincaré算子的特征值，文獻[14]通過坐標變換將一般非周期復合域轉(zhuǎn)換為規(guī)則域，并對彈性特征值問題進行了多尺度漸近分析和計算等等.

本文研究了擬周期性孔洞材料模型的橢圓型特征值問題.我們先進行理論分析，然后給出數(shù)值算例，驗證了二階雙尺度方法在解決曲線坐標下多孔材料模型特征值問題的有效性. 除非特別指出，本文中均使用Einstein求和約定，即重復指標表示求和.

2 特征值問題及坐標變換

為簡單起見，我們僅考慮二維問題. 如圖1所示，區(qū)域Ωε在直角坐標系下是擬周期區(qū)域，通過選擇適當變換ξ=L(x)，可使得擬周期區(qū)域Ωε中的每個點x唯一地映射為周期區(qū)域Θε中的點ξ. 于是可以在周期區(qū)域Θε中用漸近展開方法對特征值和特征函數(shù)進行展開，并在Θε對應的單胞區(qū)域Q*上定義單胞函數(shù)，得到問題的均勻化解、均勻化系數(shù)和多尺度逼近解，最后通過逆變換x=L-1(ξ)得到原問題的解.

圖1 坐標變換示意圖Fig.1 Coordinate transformation diagram

對于變換后的周期孔洞區(qū)域，我們做出如下假設：假設均勻化凸區(qū)域Θ為一個有界連通開子集，有Lipschitz邊界?Θ. 定義Θ上孔洞區(qū)域為Tε，Θε=Θ-Tε稱為孔洞區(qū)域. 單胞Q=[0,1]N，單胞Q*=Q-T，其中孔洞T?Q為有界開集.稱Q*為周期孔洞區(qū)域Θε對應的參考單胞.

假設孔洞邊界滿足Neumann條件，考慮(uε,λε)對應的橢圓型特征值問題為

(1)

α1,α2>0,l=(li)1≤i≤N∈RN.

由偏導數(shù)運算

可將原方程改寫為

(2)

3 二階雙尺度分析

由二階雙尺度漸進展開方法，我們可以得到(uε,λε)的漸近展開形式

uε(ξ)=u0(ξ,η)+εu1(ξ,η)+ε2u2(ξ,η)+O(ε3),

λε=λ0+ελ1+ε2λ2+O(ε3).

將偏導數(shù)運算

以及漸近展開式代入式(2)，比較方程兩端ε各冪次系數(shù)，得到

(3)

(4)

(5)

其中

根據(jù)式(3)可知，u0獨立于η.所以

u0(ξ,η)=u0(ξ)

(6)

根據(jù)式(6)，可以將式(4)化為

進而可以得到一階校正項表達式

(7)

Nα1(ξ,η)滿足單胞問題

對于式(5)，對等式兩端在Q*上做積分平均，結(jié)合式(7)，可以得到均勻化解(u0,λ0)對應均勻化方程

(8)

其中均勻化系數(shù)

將式(6)～(8)代入式(5)，整理和計算得到二階校正項表達式

Nα1α2(ξ,η)滿足單胞問題

下面考慮特征值的近似解. 基于Kesavan[10]提出的“校正方程”的思想，我們首先令wε(ξ)為以下問題的解：

相應的變分形式為

通過上述連續(xù)枯水年的供水風險分析，本論證單元在1999—2002年連續(xù)枯水年條件下，可以保證規(guī)劃水平年2015年、2020年的供水需求。論證區(qū)域2015年、2020規(guī)劃水平年的開采量，分別占多年平均地下水補給量99.3%和98.9%，雖然典型年和連續(xù)枯水年須分別動用地下水儲存量的3.4%和4.4%，但動用的地下水儲存量在豐水年、平水年可以得到回補，因此供水水源可靠。

根據(jù)

我們得到校正方程

由二階雙尺度方法，可將wε(ξ)漸近展開為

wε(ξ)=

u0(ξ)+εw1(ξ,η)+ε2w2(ξ,η)+O(ε3).

將uε(ξ),wε(ξ),λε展開式代入校正方程，比較方程兩端ε各冪次系數(shù)，可以得到

根據(jù)均勻化理論[12]，假設uε(ξ),wε(ξ)充分靠近.則

4 有限元計算

其中

當變換為線性變換時，gij全為常數(shù)，且J為正，hi=0. 此外，所有單胞函數(shù)獨立于ξ，于是一、二階逼近解為

λε,1=λ0+ελ1,λε,2=λε,1+ε2λ2.

這里Nα1(η)滿足單胞問題

Nα1α2(η)滿足單胞問題

因均勻化解滿足

于是各單胞問題以及均勻化問題的變分形式為

?φ∈H1(Q*;?Q)，

其中H1(Q*,?Q)表示空間H1(Q*,?Q)={φ|φ∈H1(Q*),φ=0 on ?Q}.

二階雙尺度有限元法的一般計算流程如下:

(1) 設定材料參數(shù)，初始化單胞區(qū)域網(wǎng)格Q*與均勻化網(wǎng)格Θ，計算函數(shù)gijkl,J,hijl.

(3) 利用求得的均勻化解與一階、二階單胞函數(shù)組裝得到一階和二階逼近解(uε,1,λε,1),(uε,2,λε,2).

(4) 將特征函數(shù)u0,uε,1,uε,2映射到原宏觀區(qū)域Ωε.

本節(jié)最后給出誤差分析方法. 首先，引入空間L2(Ωε)中內(nèi)積與范數(shù)

由于任何特征函數(shù)乘以一個非零實數(shù)也是原問題的特征函數(shù)，因而求得的特征函數(shù)不能直接進行比較. 為了進行誤差分析，我們利用最小二乘思想選擇使得兩個特征函數(shù)在L2范數(shù)中誤差最小的比例因子，且考慮縮放漸近解.也就是說，選擇比例因子使得兩個特征函數(shù)在L2范數(shù)中有最小的誤差. 為此，對于單特征值對應的特征函數(shù)，我們定義其相對誤差為

由此得到比例因子

其中

注意上式中也沒有對n求和.

5 數(shù)值算例

考慮具有復合結(jié)構(gòu)的孔洞材料，兩種材料都是均勻且各向同性的. 通過變換

x1=ξ1+ωξ2,x2=ξ2,

容易得到原宏觀區(qū)域Ωε、周期區(qū)域Θε和單胞區(qū)域Q*如圖2圖3所示.

圖2 ε=1/8時原宏觀區(qū)域Ωε和變換后區(qū)域Θε

圖3 單胞區(qū)域Q*

根據(jù)分析得

g11=1+ω2,g12=g21=-ω,g22=1,

J=1,h1=h2=0,

變換后問題變?yōu)?/p>

假設材料系數(shù)為α1=1,α2=0.01，變換系數(shù)為ω=0.25. 根據(jù)上式，應用有限元方法可得到Nα1(η)，可以計算出均勻化系數(shù)

進一步計算可得到Nα1α2(η)，從而得到一階和二階逼近解. 最終通過逆變換得到原問題的解. 由于原問題精確解難以求出，這里我們將在原宏觀區(qū)域細網(wǎng)格下計算得到的解作為精確解與各階逼近解進行比較.

表1顯示了ε=1/8和ε=1/16時的網(wǎng)格信息.隨ε減小，原宏觀細網(wǎng)格規(guī)模增大，網(wǎng)格尺寸減小，原宏觀細網(wǎng)格有限元計算量增加. 使用雙尺度有限元方法，單胞網(wǎng)格和均勻化網(wǎng)格的單元數(shù)與結(jié)點數(shù)保持不變，計算時間固定.當區(qū)域周期ε=1/16時，單胞個數(shù)增多，雙尺度有限元法更具計算優(yōu)勢.

表1 網(wǎng)格信息

下面首先對特征值進行分析.ε=1/8的前10個特征值和相對誤差如表2所示. 從表中看出，均勻化和一階校正的結(jié)果之間相差并不大，與實際結(jié)果之間則差距較大.經(jīng)過二階校正，所得到的結(jié)果明顯好于一階，因而進行二階校正很有必要.

表2 ε=1/8時前10個特征值和相對誤差

下面我們對特征函數(shù)進行分析. 圖4和圖5分別顯示了最小特征值(單特征值)和第二個特征值(2重特征值)對應細網(wǎng)格下的特征函數(shù)以及各階逼近特征函數(shù)的圖像. 從圖中容易看出，均勻化特征函數(shù)圖像光滑且與精確解相差較大，一階逼近解圖像出現(xiàn)局部振蕩，但和精確解仍有很大差別，而經(jīng)過二階校正，二階逼近解圖像已經(jīng)和精確解很接近了.

圖4 ε=1/8時第一個特征值對應的特征函數(shù)

圖5 ε=1/8時第二個特征值對應的特征函數(shù)

表3給出了L2范數(shù)下前10個特征值對應的特征函數(shù)的相對誤差.能夠看出，均勻化結(jié)果與精確解誤差較大，通過一階校正能夠有效減小誤差，但經(jīng)過一階校正后仍然不能得到滿意的估計結(jié)果，二階校正的精度遠遠大于一階. 此外，雖然進行二階校正能夠有效減小誤差，但隨著n增加特征值增大，所對應特征函數(shù)的相對誤差也越大.

表3 L2范數(shù)下特征函數(shù)的相對誤差

為考慮特征函數(shù)的收斂情況，將ε減小一半，對ε=1/16時的情況進行分析. 圖6顯示了在不同周期下第10個特征值對應的特征函數(shù)圖像.容易看出，當ε減小時，精確解和二階逼近解的特征函數(shù)圖像更相近. 從計算效果來看，當ε=1/16時，進行相同的有限元計算無論是單特征值或是重特征值所對應的特征函數(shù)的相對誤差均變小了.

圖6 第十個特征值對應的特征函數(shù)

6 結(jié)論與展望

本文通過坐標變換方法發(fā)展了多孔材料結(jié)構(gòu)模型特征值問題的二階雙尺度漸近分析方法，給出了特征函數(shù)與特征值的漸近展開表示式，提出了在曲線坐標系下的有限元算法.數(shù)值結(jié)果驗證了二階雙尺度方法解決復雜多孔結(jié)構(gòu)模型特征值問題的有效性.我們的研究表明：

(1) 均勻化解僅反映問題解的宏觀性質(zhì)，有必要借助一階與二階校正反映解的局部性質(zhì). 但鑒于一階逼近解遠不如二階逼近解來得精確，為得到更好的結(jié)果，必須添加二階校正項捕捉材料在區(qū)域內(nèi)變化的局部振蕩行為，所以將逼近解擴展至二階是合理且必要的；

(2) 選取適當?shù)摩?，對于特征值，二階校正效果明顯優(yōu)于一階校正；對于特征函數(shù)，經(jīng)過二階校正，近似解與精確解已經(jīng)非常接近；

(3) 雙尺度有限元方法所使用的單胞網(wǎng)格與均勻化網(wǎng)格總是相對粗糙的，計算時間固定，當ε越趨近于0時，雙尺度有限元方法與傳統(tǒng)方法相比更具計算優(yōu)勢；

(4) 數(shù)值實驗部分考慮了前10個特征值以及單特征值對應的特征函數(shù)的相對誤差，當特征值越來越大時，使用二階雙尺度有限元方法所得結(jié)果的誤差也越來越大，因而研究特征值和特征函數(shù)的漸近性態(tài)，對尋找適當?shù)母唠A解的校正具有重要意義.