無(wú)摩擦彈性接觸問(wèn)題的自適應(yīng)交替方向乘子法

袁 欣, 張守貴

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

0 引 言

彈性接觸問(wèn)題廣泛存在于機(jī)械工程、土木工程等領(lǐng)域,比如堤壩的結(jié)合[1-3]等．彈性接觸問(wèn)題最大的特點(diǎn)就是具有很強(qiáng)的非線性性,這就使得理論分析和數(shù)值計(jì)算難以進(jìn)行,一般只能通過(guò)一些適當(dāng)?shù)臄?shù)值方法求得近似解,如常用的Uzawa算法[4-5]、交替方向乘子法[6-8]、 增廣Lagrange法[9-10]和投影算法[11-12]等．這些方法利用變分不等式理論,將非線性彈性接觸問(wèn)題轉(zhuǎn)化為一個(gè)對(duì)應(yīng)的線性變分形式進(jìn)行迭代求解．本文研究的是一個(gè)由區(qū)域Ω代表的線性彈性體和一個(gè)剛性基礎(chǔ)支撐的彈性接觸問(wèn)題,其中Ω∈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))為位移場(chǎng)．根據(jù)Hooke定律可知,應(yīng)力張量和位移場(chǎng)的線性關(guān)系為

σ(u)=C(u),

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

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

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

考慮彈性接觸問(wèn)題[7-8],即尋找位移場(chǎng)u滿(mǎn)足以下條件:

(1)

1 增廣Lagrange乘子法

為了方便后面敘述,先引入以下記號(hào):

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

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

雙線性形式為

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

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

(2)

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

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

(3)

其中雙線性形式在V×V上滿(mǎn)足以下條件:

α的選取不依賴(lài)于u的選取,且均大于零．

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

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

(4)

un-p=0．

(5)

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

(6)

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

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

(7)

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

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

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

第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步 計(jì)算輔助變量pk+1∈L2(ΓC):

(11)

第3步 計(jì)算位移uk+1∈V:

(12)

第4步 計(jì)算Lagrange乘子:

(13)

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

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

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

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

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

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

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

(14)

第3步 計(jì)算位移uk+1∈V:

(15)

第4步 計(jì)算Lagrange乘子:

(16)

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

(17)

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

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

定義以下符號(hào):

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

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

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

其中{u,λ}是問(wèn)題(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)

因?yàn)樵谶吔绂上有λ≥0,un-g≤0,(un-g,λ)ΓC=0成立,由此可得到以下兩個(gè)不等式:

(22)

(23)

由式(22)、(23)可得

(24)

由式(14)有

代入式(24)可以得到

化簡(jiǎn)可得

變換化簡(jiǎn)得到

(25)

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

即

整理等式兩邊后得到

兩邊同時(shí)取范數(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)

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

(29)

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

(30)

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

(31)

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

所以在V上,uk→u．因?yàn)樵谶吔绂上有p=un,則

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

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

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

其中

即序列ρk滿(mǎn)足

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

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

4 算 例 分 析

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

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

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

圖1 形變后的彈性體及障礙函數(shù)(算例1)圖2 g=0.01時(shí)的算法迭代次數(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時(shí)3種算法的CPU運(yùn)行時(shí)間

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

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

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

5 結(jié) 論

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