楊小娟，田萬鵬，熊勇剛*

(1.湖南工業(yè)大學(xué)理學(xué)院，中國 株洲 412007；2.湖南都市職業(yè)學(xué)院機電工程系，中國 長沙 410137)

基于MLNNI法的正交各向異性復(fù)合材料參數(shù)識別的逆算法

楊小娟1，田萬鵬2，熊勇剛1*

(1.湖南工業(yè)大學(xué)理學(xué)院，中國 株洲 412007；2.湖南都市職業(yè)學(xué)院機電工程系，中國 長沙 410137)

為有效確定二維各向異性材料的材料參數(shù),提出一種基于無網(wǎng)格局部自然鄰近插值法(MLNNI)的識別方法.該方法只需在目標域上構(gòu)建節(jié)點數(shù)組.而其逆問題即求解以模擬數(shù)值和實測數(shù)據(jù)之間偏差為目標函數(shù)的最小值,并采用復(fù)變函數(shù)微分法(CVDM)計算用于獲取新的參數(shù)的靈敏度系數(shù).區(qū)別于有限差分法,復(fù)變函數(shù)微分法對步長大小不敏感,而且若步長足夠小,靈敏度系數(shù)的精度可以非常精確.數(shù)值算例表明所提出的方法有效.

無網(wǎng)格法;參數(shù)識別;各向異性;復(fù)變函數(shù)微分法

由于強度和剛度的特性，正交各向異性復(fù)合材料在工程結(jié)構(gòu)中得到了廣泛應(yīng)用.對于結(jié)構(gòu)設(shè)計和維護，其材料參數(shù)非常重要[1-3].然而，從實驗室里試驗樣品得到的數(shù)據(jù)與從工程實際結(jié)構(gòu)元件上得到的數(shù)據(jù)有很大區(qū)別[2].而且，在某些情況下從一個更大的結(jié)構(gòu)中取出一個測試樣本非常困難，因為這會破壞材料的內(nèi)部一致性[4].因此，尋求材料參數(shù)的識別方法[1-6]，可以得到可靠的正交各向異性介質(zhì)的材料參數(shù).但是，到目前為止，大部分的材料參數(shù)識別方法基于網(wǎng)格數(shù)值方法得到的，如有限元法和邊界元法.

本文基于無網(wǎng)格局部自然鄰近插值法提出了一種新的逆算法來確定二維各向異性材料的彈性本構(gòu)參數(shù).廣義局部彼得洛夫-伽遼金無網(wǎng)格法是該方法的一種特殊情況.近年來無網(wǎng)格方法[7-12]的發(fā)展和應(yīng)用受到越來越多的關(guān)注.這些方法的主要特點是，變量是在一個集群的分散節(jié)點上構(gòu)建的，這不僅可以為生成網(wǎng)格減輕負擔(dān)，還可以更準確地描述不規(guī)則的幾何形狀.局部彼得洛夫-伽遼金無網(wǎng)格法由Atluri和Zhu創(chuàng)立，它是一種真正意義上的無網(wǎng)格法，因為它不需要任何元素或背景網(wǎng)格用作插值或整合.這種方法不同于其他無網(wǎng)格方法，區(qū)別在于它的控制方程滿足局部子域逐點使加權(quán)殘值為零.其弱形式主要適合于幾何形狀簡單的局部區(qū)域，因此，它被成功地廣泛應(yīng)用于工程問題，如彈性靜力學(xué)問題[7-8]、斷裂力學(xué)問題[9]和不可壓縮流體流動問題[10].但是，彼得洛夫-伽遼金法本身不滿足位移邊界條件的移動最小二乘法(MLS).為了克服這個缺點，基于自然鄰點插值(NNI)[11]和當?shù)厝跣问降木植縋etrov-Galerkin方法(稱為無網(wǎng)格局部自然鄰近插值(MLNNI)方法[12]) 被用來分析固體的壓力.這種方法顯示了巨大的優(yōu)勢，因為它不僅可以保持彼得洛夫-伽遼金無網(wǎng)格法的突出特點，還具有易于收斂的自然鄰近插值的本質(zhì)邊界條件的優(yōu)點[11-12].

對于本文所考慮的逆問題，采用Levenberg-Marquardt(LM)方法來迭代求解目標函數(shù)的最小值，并確定二維正交各向異性介質(zhì)中的未知材料參數(shù).作為一種基于梯度的方法，LM方法需要計算靈敏度系數(shù)，例如，參數(shù)的目標函數(shù)的導(dǎo)數(shù).雖然有限差分法(FDM)是計算靈敏度系數(shù)的最簡單方法，但它的擾動步長的選擇必須有實驗依據(jù).目前，復(fù)變函數(shù)微分法(CVDM)，作為一種數(shù)值微分法吸引了更多研究者[6,13-14]的關(guān)注.該復(fù)變函數(shù)微分法的主要優(yōu)點在于，它避免了相近數(shù)值之間的減法，因此不涉及與小步長相關(guān)聯(lián)的精度問題.本文中，用復(fù)變函數(shù)微分法結(jié)合上述LM方法精確計算靈敏度系數(shù).兩個數(shù)值例子驗證了該方法的有效性和準確性.

1 無網(wǎng)格局部自然鄰近插值方法

如上所述，無網(wǎng)格局部自然鄰近插值是用來分析工程問題.Cai和Zhu[12]首次提出了線性彈性分析方法，以減少計算成本，簡化了基本邊界條件的施加.

1.1 自然鄰近插值

自然鄰近插值是基于著名的Voronoi圖和Delaunay三角剖分建立起來的.在二維歐氏空間R2考慮一組節(jié)點N={x1,x2,x3,…,xM}，Voronoi圖是一個分區(qū)的空間被分解成子區(qū)域TI，在每個區(qū)域TI與節(jié)點xI相關(guān)聯(lián).用數(shù)學(xué)術(shù)語表示，Voronoi多邊形TI被定義為[11]：

圖1 關(guān)于x的1階和2階泰森多邊形法細胞Fig.1 1st-order and 2nd-order Voronoi cells about x

TI={x∈R2:d(x,xI)

(1)

式中d(x,xI)是x與xI之間的距離.簡單連接的Voronoi單元有共同邊界的節(jié)點產(chǎn)生的Delaunay三角剖分.Voronoi圖和相應(yīng)的Delaunay三角剖分是雙重的.

N是x自然鄰居的數(shù)量.參考圖1，x有4個(n=4)自然節(jié)點(即節(jié)點1～4)，形狀函數(shù)φ1(x)可以表示為：φ1(x)=A1/A2,這里的A1和A2分別指四邊形abfe和abcd的面積.

關(guān)于點x的位移u(x)的一般形式可以近似寫成：

(2)

其中，uI(I=1,2,…,n)是x的n自然鄰居節(jié)點位移向量.這里需要強調(diào)的是，自然的鄰居(Sibson)形狀函數(shù)具有一些顯著的屬性，例如：積極性、插值，單位分解以及線性完整性.NNI中的更多相關(guān)細節(jié)可以參考文獻[11].

1.2 二維正交的物質(zhì)體的MLNNI公式

對于在被Γ限定的Ω的領(lǐng)域里二維線性彈性問題，控制方程和邊界條件分別是：

σij,j+bi=0, 在Ω空間里,

(3)

(4a)

(4b)

在平面應(yīng)力條件下，同類正交彈性體的本構(gòu)關(guān)系可以表示為[1-2]

σ=Dε.

(5)

其中，D=S-1是3×3彈性常數(shù)矩陣，S可以用合規(guī)管理組件寫成：

值得注意的是，只有4個五合規(guī)組件是獨立的，因為s12=s21.此外，矩形笛卡兒坐標系統(tǒng)在這里是被假定在主要材料方向的.

(6)

圖2 當?shù)囟噙呅巫佑蚝凸?jié)點I的范圍 Fig.2 Local polygonal sub-domain and boundary for node I

利用式 (6)中的散度定理，可以得到：

(7)

(8)

(9a)

(9b)

矩陣VI和BJ可表示為：

(10)

2 無網(wǎng)格局部自然鄰域插值法的逆分析公式

2.1 目標函數(shù)

這部分的目的是通過解決逆問題來確定未知的二維正交介質(zhì)的材料參數(shù).未知的材料屬性在設(shè)計向量中定義為：p=[p1,p2,…,pm]T，其中，m是未知參數(shù)的數(shù)量.在這部分，向量P的要數(shù)是正交介質(zhì)的合規(guī)組件.例如：s11,s22,s12及s66.這些未知參數(shù)可以通過下面的目標向量最小化來決定.

(11)

2.2Levenberg-Marquardt法

在這里，LM[15]法是被用來反復(fù)地最小化式(11)中的目標向量以及確定二維正交的介質(zhì)的未知參數(shù)的.LM方法繼承了高斯牛頓算法的速度優(yōu)勢以及最速下降法的穩(wěn)定性.p(k+1)=p(k)+δ(k)以一個初始可行的猜測p為起點，LM方法采用一系列的迭代步驟，根據(jù)一些特定的標準，來調(diào)整參數(shù)，直到達到集合.被調(diào)整的參數(shù)在迭代參數(shù)K的情況下以下列形式被計算[15].

(J(p(k))TJ(p(k))+μ(k)I)δ(k)=-J(p(k))Tf(p(k)),p(k+1)=p(k)+δ(k),

其中，μ(k)是正標量，也叫阻尼參數(shù).I是單位矩陣.J(p)是靈敏度系數(shù)矩陣，其被定義為

(12)

由于規(guī)定的錯誤公差η，在目前工作中，以下的這個標準被用來停止LM方法的迭代過程.‖J(p(k))f(p(k)‖<η.

2.3 約束條件

為了確保物理上的可實現(xiàn)性，必須對合規(guī)組件施加另一個不等式約束，并需要進一步進行研究.這個不等式約束可以寫成：

(13)

因此，在最小化過程中若違反了不等式約束(13)就意味著需要調(diào)整其中的一個未知參數(shù)，通過下面等式約束中的一個，調(diào)整未知參數(shù)可以很容易地被實現(xiàn)[2].

3 使用復(fù)變量微分法評價靈敏度系數(shù)矩陣

雖然在描述的式(12)中，定義了靈敏度系數(shù)矩陣J，但無法直接計算它.本文采用復(fù)變量微分方法(CVDM)來計算這個矩陣.

通過使用復(fù)變量微分法計算靈敏度衍生品，首先是由Lyness和Moler[13]提出的.在這種技術(shù)中，一個真正的函數(shù)f(x)的變量x被一個復(fù)雜的x+Ih代替.一個非常小的h(通常h=10-20)，函數(shù)f(x+Ih)可以展開成如下泰勒級數(shù)

(14)

孤立式(14)的虛數(shù)部分，得到

(15)

其中，符號“Im”代表了虛數(shù)部分.當使用有限精度算法的時候，通過確保h足夠小，延展的第三階錯誤可以被消除.因此，式(15)除以h得到：

(16)

從式(16)可以看出，使用復(fù)變量微分法的一階導(dǎo)數(shù)不受消減取消錯誤，因為它不涉及不同操作.比有限差分法(FDM)形成了一個巨大的優(yōu)勢.

(17)

其中ej表示這個單位向量第j個位置為1，其他位置為0.

4 數(shù)值算例

為了驗證該算法，在這一節(jié)中給出了兩個數(shù)值例子來識別二維正交介質(zhì)的材料參數(shù)，考慮了平面應(yīng)力狀態(tài).矩形笛卡爾坐標方向與主要材料的方向一致.在識別過程中，使用了材料參數(shù)的實際值，從無網(wǎng)格局部自然鄰域插值法中計算出來的位移是需要測量的量.材料的實際參數(shù)與其上下界在表1中列出.阻尼參數(shù)的初始值為μ(0)=10-12.

表1 材料參數(shù)的實際值和上下界

正如圖3所呈現(xiàn)的，第1個算例演示了簡支梁的受力情況.

圖3 在均布荷載條件下的正交異向簡支梁Fig.3 An orthotropic simply beam under uniform load

簡支梁上部受到均布載荷P=20 MPa.要確定的參數(shù)是4個獨立的合規(guī)組件，例如s11,s22,s12,s66.

圖4中的簡支梁被劃分為21×7個節(jié)點.在識別過程中，在水平位移上的節(jié)點1，2和3被視為測量數(shù)據(jù)，如圖4所示.材料的實際參數(shù)與預(yù)測值如表2所示.與材料的實際參數(shù)相比，識別參數(shù)是非常準確的.

圖4 正交異向簡支梁的節(jié)點排列Fig.4 Nodal arrangements for the orthotropic simply beam

實例參數(shù)/(10-10Pa-1)預(yù)測值實際值預(yù)測值/實際值確定值迭代1s110．0150．0680740．220．068074s220．250．918270．270．91827s12-0．008-0．0258680．31-0．025868s660．451．56010．291．560132s110．0150．0680740．220．068074s223．00．918273．270．91827s12-0．008-0．0258680．31-0．025868s665．51．56013．531．560153s110．30．0680744．410．068074s223．00．918273．270．91827s12-0．1-0．0258683．87-0．025868s660．451．56010．291．560154s110．30．0680744．410．068074s220．250．918270．270．91827s12-0．1-0．0258683．87-0．025868s665．51．56013．531．56015

[1] 孫秀山,黃立新,劉應(yīng)華,等. 二維正交各向異性結(jié)構(gòu)彈塑性問題的邊界元分析[J]. 復(fù)合材料學(xué)報, 2005,22(3):156-161.

[2] HUANG L X, SUN X S, LIU Y H,etal. Parameter identification for two-dimensional orthotropic material bodies by the boundary element method[J]. Engi Anal Bound Elements, 2004,28(3):109-121.

[3] 黃立新，向志海，孫秀山，等.正交各向異性復(fù)合材料孔板性能參數(shù)識別測點的最優(yōu)布置[J]. 工程力學(xué), 2006,38(10):201-209.

[4] MEUWISSEN M H, OMENS C W, BAAIJENS F P,etal. Determination of the elasto-plastic properties of aluminium using a mixed numerical-experimental method[J]. J Mater Proc Tech, 1998,75(8):204-211.

[5] LIU G R, HAN X. Computational inverse techniques for nondestructive evaluation [M]. Florida: CRC Press, 2003.

[6] GAO X W, HE M C. A new inverse analysis approach for multi-region heat conduction BEM using complex-variable-differentiation method[J]. Engi Anal Bound Elements, 2005,29(10):788-795.

[7] HU D A, LONG S Y, LIU K Y,etal. A modified meshless local Petrov-Galerkin method to elasticity problems in computer modeling and simulation[J]. Engi Anal Bound Elements, 2006,30(10):399-404.

[8] ARBDOLLAHIFA A, NAMI M R, SHAFIEI A R. A new MLPG method for elastostatic problems [J]. Engi Anal Bound Elements, 2012,36(12):451-457.

[9] CHING H K, BATRA R C. Determination of crack tip fields in linear elastostatics by the meshless local Petrov-Galerkin (MLPG) method[J]. Comp Model Engi Sci, 2001,2(5):273-289.

[10] WU X H, TAO W Q, SHEN S P,etal. A stabilized MLPG method for steady state incompressible fluid flow simulation[J]. J Comput Phy, 2010,229(12):8564-8577.

[11] SUKUMAR N, MORAN B, BELYTSCHKO T. The natural element method in solid mechanics[J]. Int J Numer Meth Engi, 1998,43(11):839-887.

[12] CAI Y C, ZHU H H. A meshless local natural neighbour interpolation method for stress analysis of solids [J]. Engi Anal Bound Elements, 2004,28(8):607-613.

[13] LYNESS J N, MOLER C B. Numerical differentiation of analytic functions [J]. Siam J Numer Anal, 1967,4(10):202-210.

[14] VATSA V N. Computation of sensitivity derivatives of navier-stokes equations using complex variables [J]. Adv Engi Softw, 2000,31(6):655-659.

[15] 袁亞湘,孫文瑜. 最優(yōu)化理論與方法[M]. 北京：科學(xué)出版社,2001.

(編輯 HWJ)

An Inverse Approach to Identify Material Properties of an Orthotropic Medium by the MLNNI Method

YANGXiao-juan1,TIANWang-peng2,XIONGYong-gang1*

(1.College of Science, Hunan University of Technology, Zhuzhou 412007, China; 2. Department of Mechanical and Electrical Engineering, Hunan Urban Professional College, Changsha 410137, China)

To determine material parameters of two-dimensional orthotropic materials, a new identification approach is proposed in this work based on the meshless local natural neighbor interpolation (MLNNI) method. It is only necessary to construct an array of nodes in the targeted domain. The identification inverse problem is formulated as the minimization of an objective function representing the difference between numerical simulation displacements and measured data. Sensitivity coefficients used to obtain parameter updates are calculated by the complex variable differentiation method (CVDM). Unlike the finite difference method, CVDM has the advantage of step size insensitivity and sufficiently small steps. The accuracy of the sensitivity coefficients can approach the computer precision. Based on the investigation of numerical example, high accuracy results have been obtained, which demonstrates the potential of our proposed approach.

meshless method; parameter identification; orthotropic; complex variable differentiation method

10.7612/j.issn.1000-2537.2017.04.011

2016-07-27

湖南省自然科學(xué)基金資助項目(2017JJ2065)；教育部人文社會科學(xué)研究青年基金資助項目(10YJC630338)

TU313.3

A

1000-2537(2017)04-0062-06

*通訊作者,E-mail：xygyxj@163.com