求解固體力學(xué)反問題的廣義邊界控制方法

王明清，劉桂榮

（1.浪潮集團高效能服務(wù)器和存儲技術(shù)國家重點實驗室，山東 濟南250101；2.太原理工大學(xué) 應(yīng)用力學(xué)與生物醫(yī)學(xué)工程研究所，山西 太原030024；

3.美國辛辛那提大學(xué) 航空航天學(xué)院，美國 俄亥俄州45221）

0 引 言

對于固體力學(xué)而言，通常需要控制方程以及所有邊上的邊界條件，進而求出問題域位移，應(yīng)力以及應(yīng)變等未知量.一般來說，只要物理問題是適定的，那么問題的解是存在且唯一的.然而，在實際的工程應(yīng)用以及科學(xué)研究中，很大一部分問題是不具有這種適定性的［1］.對于上述固體力學(xué)問題，當(dāng)所有邊界上的信息都已知，其適定性才可得到保障.但有些情況下，邊界條件是不可測的［1］.例如，在懸臂梁模型中，固定端位移的測量是困難的，同時自由端位移及應(yīng)力的測量是比較容易的.因此，常規(guī)求解正問題的方法是不能得到穩(wěn)定精確的數(shù)值解的.在同一邊界上同時給出了兩類邊界條件，這是一類典型的柯西問題.這類柯西反問題是不適定的，即測量數(shù)據(jù)微小的擾動將導(dǎo)致結(jié)果產(chǎn)生巨大誤差［1］.Lars Eldén［2］等人已經(jīng)給出了關(guān)于橢圓型偏微分方程柯西問題的穩(wěn)定性分析.其他微分方程的柯西問題可參見文獻［3-8］.

基于有限元方法求解柯西反問題的精確有效且穩(wěn)定的數(shù)值方法，首次引入了廣義邊界控制的思想，將問題域未知邊界上位移的求解轉(zhuǎn)化成求解包含一組參數(shù)化后未知系數(shù)的線性代數(shù)方程組的直接方法，不同于其他方法，無需多次迭代，格式更加簡潔.Tikhonov正則化有效地克服了隨機擾動帶來的不適定性.

本文介紹了固體力學(xué)柯西反問題，提出了有限元方法求解固體力學(xué)反問題的廣義邊界控制格式，并給出數(shù)值例子，驗證了算法的有效性和穩(wěn)定性.

1 問題的提出

1.1 正問題

如圖1所示的二維懸臂梁，其區(qū)域為Ω，邊界為Γ，其中Γ1∪Γ2∪Γ3∪Γ4=Γ且Γi∩Γj=φ，（i≠j）.懸臂梁材料為各向同性的彈性材料，施加在懸臂梁右邊界Γ4上的牽引力記為fb=［fx，fy］T，而在其余的邊界上沒有牽引力.其中，fx為x方向上的力，fy為y方向上的力.Γ1上的位移是給定的，可以記為UΓ1=［uΓ1，vΓ1］T.

相容方程為

本構(gòu)方程（胡克定律）為

正問題的控制方程可由動力學(xué)平衡方程給出，其矩陣格式表述為

其邊界條件為

圖1 右端施加牽引力的懸臂梁模型 Fig.1 A cantilever beam subjected to a traction on the right edge

1.2 柯西反問題

懸臂梁結(jié)構(gòu)是工程系統(tǒng)中最常用的部件之一.然而，直接表述其左邊界Γ1上的位移通常是非常困難的，多數(shù)情況下以UΓ1=0.因此，可得到前面定義的過于理想化的數(shù)學(xué)正問題.這種過度的簡化往往會造成應(yīng)力分析中巨大誤差.為了使應(yīng)力分析更加精確，不得不測量出邊界Γ1的實際位移，使正問題具有更加可靠的解.但直接測量出邊界Γ1上的位移通常是比較困難的，而測量邊界Γ4是非常方便的.換言之，Γ1上不給任何條件，Γ4上既給出Dirichlet邊界條件又給出Neumann邊界條件，這是一個柯西反問題

這類柯西問題是不適定的9：測量數(shù)據(jù)微小的擾動將導(dǎo)致結(jié)果產(chǎn)生巨大誤差.

2 基于有限元方法的廣義邊界控制

鑒于有限元方法是求解彈性懸臂梁正問題（3）～（6）有效的方法，本文求解反問題的過程是基于有限元方法提出的.

首先將問題域離散成一組單元與節(jié)點，然后構(gòu)造形函數(shù)，通過Galerkin弱形式構(gòu)造單元剛度矩陣，組裝整體剛度矩陣，最后施加邊界條件，得到以下有限元方程

式中：K為整體剛度矩陣；U以及F分別為Ω∪Γ上所有節(jié)點的位移向量及節(jié)點力向量.

將Γ1上所有的節(jié)點從1到l編號，其它節(jié)點編號為l+1到n.重新排列矩陣或向量中的元素，等式（10）寫成分塊形式如下

式中：I為單位矩陣；矩陣O的所有元素都為零，并且

式中：ui和vi分別為水平方向及豎直方向的節(jié)點位移；l為邊界Γ1上的節(jié)點數(shù)目；（n-l）為懸臂梁有限元模型其它部分的節(jié)點數(shù)目.

通過邊界條件（4），可得F1=U1.對于等式（7）～（9）定義的反問題，邊界Γ1上的節(jié)點位移UΓ1是未知的.首先將邊界Γ1上的條件參數(shù)化：

式中：φi（x）是基函數(shù)，這里選取多項式為基函數(shù).等式（11）中的F1可以表示為

c

i及di（i=1，2，…，k）為需要求解的系數(shù)，總共有2k個未知量.

施加邊界條件后的整體剛度矩陣K是可逆的，將K-1記為

令B1=［I，-P-1Q］T，B2=［O，P-1］T，由式（10），（14）可知，

在式（15）中，B1，B2及F2是已知的，但F1是未知的.由式（12）和（13）可知，可表示為

將式（16）代入式（15）可得

為簡化描述，標記N=B1M以及η=B2F2.式（17）可重新寫為

在式（15）中，節(jié)點位移向量U是未知的.但邊界Γ4上的節(jié)點位移向量是已知的，經(jīng)過適當(dāng)?shù)呐判蚝罂梢詫懗桑?/p>

抽取等式（18）中矩陣N以及向量η中與邊界Γ4上節(jié)點編號一致的行，形成新的矩陣和向量分別為

邊界Γ4上的位移可以表示為

3 正則化方法

由于反問題的不適定性，直接運用傳統(tǒng)的高斯消去法或最小二乘法求解等式（20）是十分不穩(wěn)定的［10］，需要將不適定的問題正則化.

Tikhonov和Arsen（1986）提出的Tikhonov正則化便是最常用的方法之一.這種簡潔的非迭代方法可求得到柯西反問題的穩(wěn)定解.另外一種有效的方法是奇異值分解方法（SVD）［11］，然而對于大型問題SVD方法往往需要很大的花費［10］.

本文將采用Tikhonov正則化方法.線性代數(shù)方程組（20）的Tikhonov正則化解為

其中，fσ表示Tikhonov泛函，定義為

其中，‖·‖代表歐幾里得范數(shù)且σ為正則化參數(shù).

為得到Tikhonov正則化解Θσ，需要求出▽fσ（Θ）=0.等式（21）的解便是正則化方程的解：

當(dāng)σ=0時，便是常見的最小二乘解或高斯消元所得到的解，這樣的解往往是不穩(wěn)定的.通過增加一個正則化參數(shù)σ，可以將范數(shù)‖Θ‖控制在一個合理的范圍內(nèi).

正則化參數(shù)σ的選取仍然是一個重要的研究課題.目前，有兩個比較常用的選取方法，即GCV方法以及L曲線法.在本文的討論中主要采用L曲線法，基于一種無擾動準則確定一個合適的σ值.Hansen和O'Leary［12-13］是最早應(yīng)用L曲線法研究不同參數(shù)對于正則化系統(tǒng)的影響，L曲線法選取正則化參數(shù)的過程可以概括為：

定義一條曲線

這條曲線被稱為L曲線，合適的正則化參數(shù)σ往往是曲線“拐角”處的點所對應(yīng)的值.

本文應(yīng)用了Hansen［14］提供的Matlab代碼求解離散的不適定方程組（20）.將L曲線法選取的正則化參數(shù)記為σ＊，方程組（20）的正則化解記為Θσ＊.將Θσ＊代入等式（16）可得

最終反問題（7）～（9）轉(zhuǎn)化成了正問題（3）～（6）.其解為

4 數(shù)值驗證

為直觀驗證算法的有效性以及穩(wěn)定性，定義相對均方根誤差為

為簡化討論，考慮問題域Ω=｛（x，y）∶0≤x≤L，-5≤y≤5｝，其中L為給定常數(shù).如圖2所示，稱邊界Γ1的兩個端點分別為“up”以及“l(fā)o”節(jié)點.

圖2 懸臂梁模型問題域 Fig.2 The solution domain of the problem

例 假設(shè)材料的楊氏模量為E=3×107且泊松比為ν=0.3，施加在右邊界Γ4上的牽引力為fb=［0，1 000］T.為了得到問題的參考解，分別給出“up”節(jié)點以及“l(fā)o”節(jié)點處的位移Uup=［0，-1×10-5］T及Ulo=［0，-8×10-2］T.其余邊界上（Γ／Γ4）的應(yīng)力設(shè)為0.在這種情況下，此問題為典型的正問題.因此，可以用有限元方法求出問題域內(nèi)任意節(jié)點的位移U，也包括邊界Γ4.考慮與以上正問題所對應(yīng)的反問題，可將邊界Γ4上的位移UΓ4作為反問題中的已知條件，并且節(jié)點位移Uup與Ulo均已知.運用本文第2節(jié)中提出的廣義邊界控制及Tikhonov正則化方法求出數(shù)值解U＊.

4.1 正則化的有效性驗證

為使測試結(jié)果更具說服力，假定固體力學(xué)模型的長度為L=20，且對邊界Γ4上的位移施加δ=5%的隨機擾動.

測試結(jié)果如預(yù)期，運用最小二乘法等非正則化的方法（令等式（23）中的σ=0）求解帶擾動的線性方程組是非常不穩(wěn)定的，結(jié)果可參見圖3.圖3中實線代表有限元方法得到的參考解，點線為最小二乘解.顯然，最小二乘法求解不適定問題是十分不穩(wěn)定的.高斯消去法等類似的非正則化方法會得出相同的結(jié)論.

圖4展示了運用Tikhonov正則化且通過L曲線法選取正則化參數(shù)σ后所得到的反問題的解.與最小二乘法等非正則化方法相比，Tikhonov正則化結(jié)果極好地吻合了參考解.因此，可斷定Tikhonov正則化克服柯西反問題不適定性的有效方法.同時，L曲線法為Tikhonov正則化選取了一個合適的參數(shù)σ.

圖3 邊界Γ1上的有限元參考解以及最小二乘解（δ=5%） Fig.3 Analytical solution with FEM and LS solution（δ=5%）

圖4 邊界Γ1上的有限元參考解以及Tikhonov正則解（δ=5%） Fig.4 Analytical solution with FEM and Tikhonov regularization solution（δ=5%）

4.2 算法的有效性及穩(wěn)定性驗證

本小節(jié)對柯西邊界上的數(shù)據(jù)施加不同水平的擾動，考察解的穩(wěn)定性.假定模型的長度為L=20，對邊界Γ4上的位移施加不同水平（直至10%）的隨機擾動，并應(yīng)用基于有限元方法廣義邊界控制求解.圖5給出了不同擾動水平下的求解結(jié)果且整個問題域上的相對均方根誤差（RSE）也在表1中給出.從圖5和表1可看出，對于不同水平的隨機擾動，數(shù)值方法是非常穩(wěn)定的.且當(dāng)擾動水平已達10%時，數(shù)值解依然很好地吻合參考解.

表1 不同擾動水平下的相對平方根誤差Tab.1 The RSE against various noise levels

求解左邊界上的位移是通過右邊界上的柯西數(shù)據(jù)控制的，數(shù)值方法的靈敏度理論上會隨著模型長度的增加而降低.接下來對算法的靈敏度進行測試.在擾動水平分別為δ=1%以及δ=5%的情況下，逐漸增加梁的長度來檢驗算法的靈敏度.表2給出了測試不同長度模型所對應(yīng)的相對均方根誤差.從表2可以看出，無論對于實驗中的哪種長寬比例，廣義邊界控制方法都是十分有效的.正如預(yù)料，隨著梁模型的增長相對誤差也是逐漸加大的.但即便長寬比高達6∶1且擾動水平為5%，數(shù)值解仍然可以較好地吻合參考解.

表2 不同長寬比下的相對均方根誤差 Tab.2 The RSE with different length-width ratio

圖5 參考解及不同擾動水平下的數(shù)值解 Fig.5 Analytical solution and numerical solutions against different noise levels

5 結(jié) 論

本文提出了一個基于有限元方法求解柯西反問題的精確有效且穩(wěn)定的數(shù)值方法.首次引入了廣義邊界控制的思想，將問題域未知邊界上位移的求解轉(zhuǎn)化成求解包含一組參數(shù)化后未知系數(shù)的線性代數(shù)方程組的直接方法.Tikhonov正則化有效地克服了隨機擾動帶來的不適定性.數(shù)值論證說明了基于有限元方法的廣義邊界控制是有效且穩(wěn)定的.由于算法具有很好的魯棒性，因此對于求解實際應(yīng)用中的反問題是非常有用的.

［1］Liu G R，Han X.Computational inverse techniques in nondestructive evalution［M］.Florida：CRC Press，2003.

［2］Eldén L，Berntsson F.A stability estimate for a Cauchy problem for an elliptic partial differential equations［J］.Inverse Probl.，2005，21（5）：1643-1653.

［3］Hon Y C，Li M.A computational method for inverse free boundary determination problem［J］.Int.J.Numer.Meth.Eng.2008，73（9）：1291-1309.

［4］Hon Y C，Li M.A discrepancy principle for the source points location in using the MFS for solving BHCP［J］.Int.J.Comput.Meth.，2009，6（2）：181-197.

［5］Hon Y C，Wei T.A fundamental solution method for inverse heat conduction problem［J］.Eng.Anal.Bound.Elem.，2004，28（5）：489-495.

［6］Ling L，Takeuchi T.Boundary control for inverse Cauchy problems of the Laplace equations［J］.CMES，2008，29：45-54.

［7］Wei T，Hon Y C，Ling L.Method of fundamental solutions with regularization techniques for Cauchy problems of elliptic operators［J］.Eng.Anal.Bound.Elem.，2007，31（4）：373-385.

［8］Zhang H W，Wei T.Two iterative methods for a Cauchy problem of the elliptic equation with variable coefficients in a strip region［J］.Numer Algorithms，2014，65（4）：875-892.

［9］Schaback R.Solving the Laplace equation by meshless collocation using harmonic kernles［J］.Adv.Comput.Math.，2009，31（4）：457-470.

［10］Lawson C L，Hanson R J.Solving least squares problems［M］.Philadephia：SIAM，1995.

［11］Kalman D.Singular valuable decomposition：the SVD of a matrix［J］.The College Mathematics Journal，1996，27：2-23.

［12］Hansen P C.The L-curve and its use in the numerical treatment of inverse problems［J］.Computational Inverse Problems in Electrocardiology，2001：119-142.

［13］Hansen P C，O′Leary D P.The use of the L-curve in the regularization of discrete ill-posed problems［J］.SIAM J.Sci.Comput.，1993，14（6）：1487-1503.

［14］Hansen P C.Regularization tools：a matlab package for analysis and solution of discrete ill-posed problems［J］.Numer Algorithms，1994，6（1）：1-35.