葉舒愉， 馬 強(qiáng)， 胡 兵

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

1 引 言

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

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

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

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

2 特征值問題及坐標(biāo)變換

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

圖1 坐標(biāo)變換示意圖Fig.1 Coordinate transformation diagram

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

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

(1)

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

由偏導(dǎo)數(shù)運算

可將原方程改寫為

(2)

3 二階雙尺度分析

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

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

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

將偏導(dǎo)數(shù)運算

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

(3)

(4)

(5)

其中

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

u0(ξ,η)=u0(ξ)

(6)

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

進(jìn)而可以得到一階校正項表達(dá)式

(7)

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

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

(8)

其中均勻化系數(shù)

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

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

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

相應(yīng)的變分形式為

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

根據(jù)

我們得到校正方程

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

wε(ξ)=

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

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

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

4 有限元計算

其中

當(dāng)變換為線性變換時，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è)定材料參數(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ù)不能直接進(jìn)行比較. 為了進(jìn)行誤差分析，我們利用最小二乘思想選擇使得兩個特征函數(shù)在L2范數(shù)中誤差最小的比例因子，且考慮縮放漸近解.也就是說，選擇比例因子使得兩個特征函數(shù)在L2范數(shù)中有最小的誤差. 為此，對于單特征值對應(yīng)的特征函數(shù)，我們定義其相對誤差為

由此得到比例因子

其中

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

5 數(shù)值算例

考慮具有復(fù)合結(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è)材料系數(shù)為α1=1,α2=0.01，變換系數(shù)為ω=0.25. 根據(jù)上式，應(yīng)用有限元方法可得到Nα1(η)，可以計算出均勻化系數(shù)

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

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

表1 網(wǎng)格信息

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

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

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

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

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

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

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

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

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

6 結(jié)論與展望

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

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

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

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

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