無摩擦彈性接觸問題的自適應交替方向乘子法

袁 欣, 張守貴

(重慶師范大學 數學科學學院, 重慶 401331)

0 引 言

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

σ(u)=C(u),

其中C=(Cijkl)是(四階)彈性模張量,它是對稱正定的．設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},

雙線性形式為

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

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

(2)

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

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

(3)

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

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

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

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

(4)

un-p=0．

(5)

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

(6)

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

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

(7)

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

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

由文獻[13-14]知,可以通過改變輔助變量和位移場的計算順序對算法進行優(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步 給定某種判定條件,若滿足則停止迭代得到數值解uk+1; 否則,置k=k+1,返回第2 步．

2 自適應交替方向乘子法

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

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

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

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

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

(14)

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

(15)

第4步 計算Lagrange乘子:

(16)

第5步 選取罰參數ρk+1,使得

(17)

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

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

定義以下符號:

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

根據算法的原理和問題的性質,可得算法2的收斂結果．

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

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

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

兩式相減可得

即

(18)

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

(19)

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

(20)

而且由式(16)有

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

(21)

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

(22)

(23)

由式(22)、(23)可得

(24)

由式(14)有

代入式(24)可以得到

化簡可得

變換化簡得到

(25)

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

即

整理等式兩邊后得到

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

即定理得證．

定理2 設{uk,pk,λk}是算法2產生的序列,則在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)

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

(29)

根據式(27)可得

(30)

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

(31)

由級數收斂的必要條件可知

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

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

對于任意k≥1,由于

所以可以得到

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

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

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

3 罰參數的自適應法則

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

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

其中

即序列ρk滿足

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

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

4 算 例 分 析

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

為了驗證算法的有效性,我們給出兩個數值算例．在數值算例中,我們取τ=2和cmax=100．對所有數值結算都采用迭代終止條件‖uk+1-uk‖≤10-5‖uk+1‖．

算例1 設Ω=(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上,彈性體Ω和剛性基礎支撐之間的障礙函數g(x)=0.01．材料常數即彈性模量和Poisson比分別為2 000和0.3．

圖1 形變后的彈性體及障礙函數(算例1)圖2 g=0.01時的算法迭代次數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所示(數據擴大為原來的10倍)．為了展示算法的性能,對3種算法(ADMM1(文獻[8]中的算法1),ADMM2(本文中的算法1),SADMM(本文中的算法2))取相同的初始罰參數和步長h=20,考察了迭代收斂速度,如圖2所示．顯然,固定罰參數的交替方向乘子法(ADMM1、ADMM2)對參數的取值非常敏感,數值過大或過小都會減緩算法的收斂速度．而SADMM的罰參數自動選取,在很大程度上解決了算法對參數取值的依賴性．

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

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

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

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

5 結 論

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