袁 欣, 張守貴

(重慶師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院, 重慶 401331)

0 引 言

彈性接觸問題廣泛存在于機(jī)械工程、土木工程等領(lǐng)域,比如堤壩的結(jié)合[1-3]等．彈性接觸問題最大的特點就是具有很強(qiáng)的非線性性,這就使得理論分析和數(shù)值計算難以進(jìn)行,一般只能通過一些適當(dāng)?shù)臄?shù)值方法求得近似解,如常用的Uzawa算法[4-5]、交替方向乘子法[6-8]、 增廣Lagrange法[9-10]和投影算法[11-12]等．這些方法利用變分不等式理論,將非線性彈性接觸問題轉(zhuǎn)化為一個對應(yīng)的線性變分形式進(jìn)行迭代求解．本文研究的是一個由區(qū)域Ω代表的線性彈性體和一個剛性基礎(chǔ)支撐的彈性接觸問題,其中Ω∈R2．區(qū)域Ω的邊界Γ(Γ=ΓD∪ΓN∪ΓC)是光滑的．假設(shè)彈性體固定在ΓD上,且mean(ΓD)>0．邊界ΓC是Γ的一部分,它也是Ω和基礎(chǔ)支撐的接觸面．ΓC和剛性基礎(chǔ)支撐之間的障礙函數(shù)用g表示．在本文中考慮很小的形變,應(yīng)變張量為(u)=(?u+?uT)/2,其中u=(u1(x),…,ud(x))為位移場．根據(jù)Hooke定律可知,應(yīng)力張量和位移場的線性關(guān)系為

σ(u)=C(u),

其中C=(Cijkl)是(四階)彈性模張量,它是對稱正定的．設(shè)n表示Ω在Γ上的單位外法向量,t表示Ω在Γ上的單位切向量,則

un=u·n,ut=u-unn,

σn(u)=(σn(u)n)·n,σt(u)=σ(u)n-σn(u)n．

考慮彈性接觸問題[7-8],即尋找位移場u滿足以下條件:

(1)

1 增廣Lagrange乘子法

為了方便后面敘述,先引入以下記號:

V={v∈H1(Ω)2,v=0, onΓD},

K={v∈V,vn-g≤0, onΓC},

雙線性形式為

根據(jù)文獻(xiàn)[7]可知,彈性接觸問題(1)可以轉(zhuǎn)化為限制性極小值問題,尋找u∈K滿足不等式:

J(u)≤J(v), ?v∈K．

(2)

由文獻(xiàn)[1,9]可引入Lagrange乘子λ∈L2(Ω),得問題(1)的變分形式:

a(u,v)=l(v)-(λ,vn),

(3)

其中雙線性形式在V×V上滿足以下條件:

α的選取不依賴于u的選取,且均大于零．

為了通過交替方向乘子法解決問題(2),需要構(gòu)造一個用增廣Lagrange函數(shù)表示的鞍點問題．根據(jù)文獻(xiàn)[13],定義集合

所以問題(2)就等價于下面的限制性極小值問題,尋找(u,p)∈V×L2(ΓC)滿足

(4)

un-p=0．

(5)

根據(jù)式(4)和(5)可以得到增廣Lagrange函數(shù):

(6)

其中ρ>0為罰參數(shù)．考慮由式(6)表示如下鞍點問題:

Lρ(u,p;μ)≤Lρ(u,p;λ)≤Lρ(v,q;λ), ?{v,q;λ}∈V×L2(ΓC)×L2(ΓC),

(7)

然后可以求解問題(2)和(7)．

由引理1和文獻(xiàn)[6-8]中的交替方向乘子法(ADMM1)可知,要想求解問題(2),只需要求解增廣Lagrange函數(shù)Lρ的鞍點{u,p;λ}．即先由給定的初始值和迭代公式計算位移,再計算輔助變量,最后更新Lagrange乘子,以此順序迭代直至達(dá)到終止條件．

由文獻(xiàn)[13-14]知,可以通過改變輔助變量和位移場的計算順序?qū)λ惴ㄟM(jìn)行優(yōu)化,并采用交替方向乘子法計算該鞍點問題．算法過程如下:

第1步 給定初始值{u0,λ0}∈V×L2(ΓC),ρ>0,置k=0．

第2步 求輔助變量pk+1∈L2(ΓC):

Lρ(uk,pk+1;λk)≤Lρ(uk,q;λk), ?q∈L2(ΓC)．

(8)

第3步 求解位移uk+1∈V:

Lρ(uk+1,pk+1;λk)≤Lρ(v,pk+1;λk), ?v∈V．

(9)

第4步 更新Lagrange乘子:

(10)

算法1 交替方向乘子法(ADMM2)

第1步 給定初始值λ0,u0,ρ>0,置k=0．

第2步 計算輔助變量pk+1∈L2(ΓC):

(11)

第3步 計算位移uk+1∈V:

(12)

第4步 計算Lagrange乘子:

(13)

第5步 給定某種判定條件,若滿足則停止迭代得到數(shù)值解uk+1; 否則,置k=k+1,返回第2 步．

2 自適應(yīng)交替方向乘子法

我們注意到算法1對任何固定的罰參數(shù)ρ>0都是無條件收斂的[14]．但是,如果罰參數(shù)過大或者過小,都會影響算法的收斂速度．為了改善算法的性能,我們對罰參數(shù)應(yīng)用了自適應(yīng)法,則得到可變罰參數(shù)ρk．接下來我們假設(shè)一個非負(fù)序列{sk},滿足

在算法1中讓罰參數(shù)在迭代過程中自動選取,可以得到如下罰參數(shù)自適應(yīng)交替方向乘子法(SADMM),由此得到變參數(shù)序列{ρk},以此序列代替固定參數(shù)ρ．

算法2 自適應(yīng)交替方向乘子法(SADMM)

第1步 給定初始值λ0,u0,ρ>0,罰參數(shù)ρ0=ρ,置k=0．

第2步 計算輔助變量pk+1∈L2(ΓC):

(14)

第3步 計算位移uk+1∈V:

(15)

第4步 計算Lagrange乘子:

(16)

第5步 選取罰參數(shù)ρk+1,使得

(17)

第6步 給定某種判定條件,若滿足則停止迭代得到數(shù)值解uk+1; 否則,置k=k+1,返回第2步．

為了證明算法2的收斂性,先給出以下引理和定理．

定義以下符號:

由算法2的第5步可知,罰參數(shù)ρk∈[ρ0/Cs,Csρ0]是有界的,設(shè)

根據(jù)算法的原理和問題的性質(zhì),可得算法2的收斂結(jié)果．

定理1 設(shè){(uk,pk),λk}是由算法2產(chǎn)生的序列,則有

其中{u,λ}是問題(3)的解．

證明令δuk=uk-u,δpk=pk-p,δλk=λk-λ,根據(jù)式(3)和(15)可得

兩式相減可得

即

(18)

在式(18)中取v=δuk+1,則

(19)

由于p=un在ΓC上成立,則代入式(19)得到

(20)

而且由式(16)有

因此由以上兩式,并根據(jù)雙線性形式a(·,·)的性質(zhì)可知

(21)

因為在邊界ΓC上有λ≥0,un-g≤0,(un-g,λ)ΓC=0成立,由此可得到以下兩個不等式:

(22)

(23)

由式(22)、(23)可得

(24)

由式(14)有

代入式(24)可以得到

化簡可得

變換化簡得到

(25)

由式(16)和p=un在ΓC上成立,可得

即

整理等式兩邊后得到

兩邊同時取范數(shù)并由式(25)可得

即定理得證．

定理2 設(shè){uk,pk,λk}是算法2產(chǎn)生的序列,則在V上uk→u,在L2(ΓC)上pk→p,在L2(ΓC)上λk?λ．證明由于sk≥0, 0<ρk+1≤(1+sk)ρk,可得

(26)

由定理1可知

從而由式(26)有

(27)

由式(21)和(27)可得

(28)

即存在一個常數(shù)C>0使得下面不等式成立:

(29)

根據(jù)式(27)可得

(30)

由式(21)、(27)和(30),結(jié)合式(21),我們可以得到

(31)

由級數(shù)收斂的必要條件可知

所以在V上,uk→u．因為在邊界ΓC上有p=un,則

所以在L2(ΓC)上有pk+1→p．

對于任意k≥1,由于

所以可以得到

結(jié)合式(21)可得序列{λk}是有界序列,可證λk?λ[14]．從而定理得證．

綜上,算法2的收斂性得證．

注1 若在算法2的第5步中令sk=0,則ρk=ρ0．因此算法2為固定罰參數(shù)情形,即為算法1,故同理可證算法1也是收斂的．

3 罰參數(shù)的自適應(yīng)法則

利用自適應(yīng)交替方向乘子法求解無摩擦單側(cè)彈性接觸問題,迭代過程中通過迭代函數(shù)自動調(diào)整罰參數(shù),用變參數(shù)ρk代替固定參數(shù)ρ[10-11,14-15],從而達(dá)到提高算法效率的目的．下面具體考慮算法2中選取罰參數(shù)ρk的自適應(yīng)法則．由定理1 的結(jié)論可知序列{uk,pk,λk}滿足以下不等式:

為了提高算法的收斂速度,利用如下平衡原理:

其中

即序列ρk滿足

序列{sk}按照如下方式得到:

其中cmax是一個正常數(shù),ck表示ρk改變的次數(shù),即

4 算 例 分 析

本文使用FreeFEM++軟件求解算例的位移場、迭代次數(shù)和CPU運行時間．該軟件可以對二維和三維偏微分問題自動建立有限元模型,我們只需要在此基礎(chǔ)上定義網(wǎng)格和問題即可進(jìn)行求解．

為了驗證算法的有效性,我們給出兩個數(shù)值算例．在數(shù)值算例中,我們?nèi)ˇ?2和cmax=100．對所有數(shù)值結(jié)算都采用迭代終止條件‖uk+1-uk‖≤10-5‖uk+1‖．

算例1 設(shè)Ω=(0,2)×(0,1)是一個矩形的彈性體,邊界由ΓD={0}×(0,1),ΓN1=(0,2)×{1},ΓN2={2}×(0,1)和ΓC=(0,2)×{0}四部分組成．彈性體在ΓD上固定不動,即u=0,在ΓN1上滿足Riemann邊界條件σ(u)n=(0,-10)T,在ΓN2上滿足σ(u)n=0．在ΓC上,彈性體Ω和剛性基礎(chǔ)支撐之間的障礙函數(shù)g(x)=0.01．材料常數(shù)即彈性模量和Poisson比分別為2 000和0.3．

圖1 形變后的彈性體及障礙函數(shù)(算例1)圖2 g=0.01時的算法迭代次數(shù)Fig. 1 The deformed configuration and the Fig. 2 The number of iterations for each obstacle function (example 1)method with g=0.01

表1 g=0.01時3種算法的CPU運行時間

用算法2求解此算例得到彈性體形變圖,如圖1所示(數(shù)據(jù)擴(kuò)大為原來的10倍)．為了展示算法的性能,對3種算法(ADMM1(文獻(xiàn)[8]中的算法1),ADMM2(本文中的算法1),SADMM(本文中的算法2))取相同的初始罰參數(shù)和步長h=20,考察了迭代收斂速度,如圖2所示．顯然,固定罰參數(shù)的交替方向乘子法(ADMM1、ADMM2)對參數(shù)的取值非常敏感,數(shù)值過大或過小都會減緩算法的收斂速度．而SADMM的罰參數(shù)自動選取,在很大程度上解決了算法對參數(shù)取值的依賴性．

取不同的初始罰參數(shù)ρ和步長h, 采用3種算法對問題進(jìn)行數(shù)值計算, 所需的CPU運行時間如表1所示．圖1和表1中數(shù)值結(jié)果均表明,自適應(yīng)交替方向乘子法收斂速度及運行速度明顯優(yōu)于固定罰參數(shù)交替方向乘子法,并且對所有的罰參數(shù)具有較好的穩(wěn)定性．

圖3 形變后的彈性體及障礙函數(shù)(算例2)圖4 各算法的迭代次數(shù)Fig. 3 The deformed configuration and the obstacle Fig. 4 The number of iterations for each method function (example 2)

表2 3種算法的CPU運行時間

對此算例應(yīng)用算法2得到形變后的彈性體,如圖3所示(數(shù)據(jù)擴(kuò)大為原來的20倍)．對3種算法(ADMM1、ADMM2、SADMM) 取不同的初始罰參數(shù)ρ和步長h=20,考察算法收斂速度,數(shù)值結(jié)果如圖4所示．從數(shù)值結(jié)果可以看出,ADMM2比ADMM1的迭代次數(shù)更少,但是還是沒有擺脫對罰參數(shù)的依賴性．而SADMM相對比較穩(wěn)定,不會因為罰參數(shù)的過大或者過小而顯著增加迭代次數(shù)．為了更全面展示SADMM的性能,分別取不同的初始罰參數(shù)ρ和步長h,考察3種算法的CPU運行時間,如表2所示．結(jié)果表明自適應(yīng)交替方向乘子法具有較好的性能．

5 結(jié) 論

本文提出了求解彈性接觸問題的自適應(yīng)交替方向乘子法,該算法不僅計算簡單,并且給出的自適應(yīng)法則利用邊界迭代函數(shù)自動選擇合適的罰參數(shù),能有效解決交替方向乘子法的迭代次數(shù)高度依賴于罰參數(shù)ρ的問題,從而顯著提高算法性能．?dāng)?shù)值結(jié)果表明了該算法的有效性．