楊小娟，田萬(wàn)鵬，熊勇剛*

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

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

楊小娟1，田萬(wàn)鵬2，熊勇剛1*

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

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

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

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

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

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

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

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

1.1 自然鄰近插值

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

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

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

(1)

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

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

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

(2)

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

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

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

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

(3)

(4a)

(4b)

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

σ=Dε.

(5)

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

值得注意的是，只有4個(gè)五合規(guī)組件是獨(dú)立的，因?yàn)閟12=s21.此外，矩形笛卡兒坐標(biāo)系統(tǒng)在這里是被假定在主要材料方向的.

(6)

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

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

(7)

(8)

(9a)

(9b)

矩陣VI和BJ可表示為：

(10)

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

2.1 目標(biāo)函數(shù)

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

(11)

2.2Levenberg-Marquardt法

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

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

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

(12)

由于規(guī)定的錯(cuò)誤公差η，在目前工作中，以下的這個(gè)標(biāo)準(zhǔn)被用來(lái)停止LM方法的迭代過(guò)程.‖J(p(k))f(p(k)‖<η.

2.3 約束條件

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

(13)

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

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

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

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

(14)

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

(15)

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

(16)

從式(16)可以看出，使用復(fù)變量微分法的一階導(dǎo)數(shù)不受消減取消錯(cuò)誤，因?yàn)樗簧婕安煌僮?比有限差分法(FDM)形成了一個(gè)巨大的優(yōu)勢(shì).

(17)

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

4 數(shù)值算例

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

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

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

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

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

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

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

實(shí)例參數(shù)/(10-10Pa-1)預(yù)測(cè)值實(shí)際值預(yù)測(cè)值/實(shí)際值確定值迭代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)彈塑性問(wèn)題的邊界元分析[J]. 復(fù)合材料學(xué)報(bào), 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ù)識(shí)別測(cè)點(diǎn)的最優(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é)基金資助項(xiàng)目(2017JJ2065)；教育部人文社會(huì)科學(xué)研究青年基金資助項(xiàng)目(10YJC630338)

TU313.3

A

1000-2537(2017)04-0062-06

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