基于MLNNI法的正交各向异性复合材料参数识别的逆算法

杨小娟，田万鹏，熊勇刚*

(1.湖南工业大学理学院，中国 株洲 412007；2.湖南都市职业学院机电工程系，中国 长沙 410137)

基于MLNNI法的正交各向异性复合材料参数识别的逆算法

杨小娟1，田万鹏2，熊勇刚1*

(1.湖南工业大学理学院，中国 株洲 412007；2.湖南都市职业学院机电工程系，中国 长沙 410137)

为有效确定二维各向异性材料的材料参数,提出一种基于无网格局部自然邻近插值法(MLNNI)的识别方法.该方法只需在目标域上构建节点数组.而其逆问题即求解以模拟数值和实测数据之间偏差为目标函数的最小值,并采用复变函数微分法(CVDM)计算用于获取新的参数的灵敏度系数.区别于有限差分法,复变函数微分法对步长大小不敏感,而且若步长足够小,灵敏度系数的精度可以非常精确.数值算例表明所提出的方法有效.

无网格法;参数识别;各向异性;复变函数微分法

由于强度和刚度的特性，正交各向异性复合材料在工程结构中得到了广泛应用.对于结构设计和维护，其材料参数非常重要[1-3].然而，从实验室里试验样品得到的数据与从工程实际结构元件上得到的数据有很大区别[2].而且，在某些情况下从一个更大的结构中取出一个测试样本非常困难，因为这会破坏材料的内部一致性[4].因此，寻求材料参数的识别方法[1-6]，可以得到可靠的正交各向异性介质的材料参数.但是，到目前为止，大部分的材料参数识别方法基于网格数值方法得到的，如有限元法和边界元法.

本文基于无网格局部自然邻近插值法提出了一种新的逆算法来确定二维各向异性材料的弹性本构参数.广义局部彼得洛夫-伽辽金无网格法是该方法的一种特殊情况.近年来无网格方法[7-12]的发展和应用受到越来越多的关注.这些方法的主要特点是，变量是在一个集群的分散节点上构建的，这不仅可以为生成网格减轻负担，还可以更准确地描述不规则的几何形状.局部彼得洛夫-伽辽金无网格法由Atluri和Zhu创立，它是一种真正意义上的无网格法，因为它不需要任何元素或背景网格用作插值或整合.这种方法不同于其他无网格方法，区别在于它的控制方程满足局部子域逐点使加权残值为零.其弱形式主要适合于几何形状简单的局部区域，因此，它被成功地广泛应用于工程问题，如弹性静力学问题[7-8]、断裂力学问题[9]和不可压缩流体流动问题[10].但是，彼得洛夫-伽辽金法本身不满足位移边界条件的移动最小二乘法(MLS).为了克服这个缺点，基于自然邻点插值(NNI)[11]和当地弱形式的局部Petrov-Galerkin方法(称为无网格局部自然邻近插值(MLNNI)方法[12]) 被用来分析固体的压力.这种方法显示了巨大的优势，因为它不仅可以保持彼得洛夫-伽辽金无网格法的突出特点，还具有易于收敛的自然邻近插值的本质边界条件的优点[11-12].

对于本文所考虑的逆问题，采用Levenberg-Marquardt(LM)方法来迭代求解目标函数的最小值，并确定二维正交各向异性介质中的未知材料参数.作为一种基于梯度的方法，LM方法需要计算灵敏度系数，例如，参数的目标函数的导数.虽然有限差分法(FDM)是计算灵敏度系数的最简单方法，但它的扰动步长的选择必须有实验依据.目前，复变函数微分法(CVDM)，作为一种数值微分法吸引了更多研究者[6,13-14]的关注.该复变函数微分法的主要优点在于，它避免了相近数值之间的减法，因此不涉及与小步长相关联的精度问题.本文中，用复变函数微分法结合上述LM方法精确计算灵敏度系数.两个数值例子验证了该方法的有效性和准确性.

1 无网格局部自然邻近插值方法

如上所述，无网格局部自然邻近插值是用来分析工程问题.Cai和Zhu[12]首次提出了线性弹性分析方法，以减少计算成本，简化了基本边界条件的施加.

1.1 自然邻近插值

自然邻近插值是基于著名的Voronoi图和Delaunay三角剖分建立起来的.在二维欧氏空间R2考虑一组节点N={x1,x2,x3,…,xM}，Voronoi图是一个分区的空间被分解成子区域TI，在每个区域TI与节点xI相关联.用数学术语表示，Voronoi多边形TI被定义为[11]：

图1 关于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单元有共同边界的节点产生的Delaunay三角剖分.Voronoi图和相应的Delaunay三角剖分是双重的.

N是x自然邻居的数量.参考图1，x有4个(n=4)自然节点(即节点1～4)，形状函数φ1(x)可以表示为：φ1(x)=A1/A2,这里的A1和A2分别指四边形abfe和abcd的面积.

关于点x的位移u(x)的一般形式可以近似写成：

(2)

其中，uI(I=1,2,…,n)是x的n自然邻居节点位移向量.这里需要强调的是，自然的邻居(Sibson)形状函数具有一些显著的属性，例如：积极性、插值，单位分解以及线性完整性.NNI中的更多相关细节可以参考文献[11].

1.2 二维正交的物质体的MLNNI公式

对于在被Γ限定的Ω的领域里二维线性弹性问题，控制方程和边界条件分别是：

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

(3)

(4a)

(4b)

在平面应力条件下，同类正交弹性体的本构关系可以表示为[1-2]

σ=Dε.

(5)

其中，D=S-1是3×3弹性常数矩阵，S可以用合规管理组件写成：

值得注意的是，只有4个五合规组件是独立的，因为s12=s21.此外，矩形笛卡儿坐标系统在这里是被假定在主要材料方向的.

(6)

图2 当地多边形子域和节点I的范围 Fig.2 Local polygonal sub-domain and boundary for node I

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

(7)

(8)

(9a)

(9b)

矩阵VI和BJ可表示为：

(10)

2 无网格局部自然邻域插值法的逆分析公式

2.1 目标函数

这部分的目的是通过解决逆问题来确定未知的二维正交介质的材料参数.未知的材料属性在设计向量中定义为：p=[p1,p2,…,pm]T，其中，m是未知参数的数量.在这部分，向量P的要数是正交介质的合规组件.例如：s11,s22,s12及s66.这些未知参数可以通过下面的目标向量最小化来决定.

(11)

2.2Levenberg-Marquardt法

在这里，LM[15]法是被用来反复地最小化式(11)中的目标向量以及确定二维正交的介质的未知参数的.LM方法继承了高斯牛顿算法的速度优势以及最速下降法的稳定性.p(k+1)=p(k)+δ(k)以一个初始可行的猜测p为起点，LM方法采用一系列的迭代步骤，根据一些特定的标准，来调整参数，直到达到集合.被调整的参数在迭代参数K的情况下以下列形式被计算[15].

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

其中，μ(k)是正标量，也叫阻尼参数.I是单位矩阵.J(p)是灵敏度系数矩阵，其被定义为

(12)

由于规定的错误公差η，在目前工作中，以下的这个标准被用来停止LM方法的迭代过程.‖J(p(k))f(p(k)‖<η.

2.3 约束条件

为了确保物理上的可实现性，必须对合规组件施加另一个不等式约束，并需要进一步进行研究.这个不等式约束可以写成：

(13)

因此，在最小化过程中若违反了不等式约束(13)就意味着需要调整其中的一个未知参数，通过下面等式约束中的一个，调整未知参数可以很容易地被实现[2].

3 使用复变量微分法评价灵敏度系数矩阵

虽然在描述的式(12)中，定义了灵敏度系数矩阵J，但无法直接计算它.本文采用复变量微分方法(CVDM)来计算这个矩阵.

通过使用复变量微分法计算灵敏度衍生品，首先是由Lyness和Moler[13]提出的.在这种技术中，一个真正的函数f(x)的变量x被一个复杂的x+Ih代替.一个非常小的h(通常h=10-20)，函数f(x+Ih)可以展开成如下泰勒级数

(14)

孤立式(14)的虚数部分，得到

(15)

其中，符号“Im”代表了虚数部分.当使用有限精度算法的时候，通过确保h足够小，延展的第三阶错误可以被消除.因此，式(15)除以h得到：

(16)

从式(16)可以看出，使用复变量微分法的一阶导数不受消减取消错误，因为它不涉及不同操作.比有限差分法(FDM)形成了一个巨大的优势.

(17)

其中ej表示这个单位向量第j个位置为1，其他位置为0.

4 数值算例

为了验证该算法，在这一节中给出了两个数值例子来识别二维正交介质的材料参数，考虑了平面应力状态.矩形笛卡尔坐标方向与主要材料的方向一致.在识别过程中，使用了材料参数的实际值，从无网格局部自然邻域插值法中计算出来的位移是需要测量的量.材料的实际参数与其上下界在表1中列出.阻尼参数的初始值为μ(0)=10-12.

表1 材料参数的实际值和上下界

正如图3所呈现的，第1个算例演示了简支梁的受力情况.

图3 在均布荷载条件下的正交异向简支梁Fig.3 An orthotropic simply beam under uniform load

简支梁上部受到均布载荷P=20 MPa.要确定的参数是4个独立的合规组件，例如s11,s22,s12,s66.

图4中的简支梁被划分为21×7个节点.在识别过程中，在水平位移上的节点1，2和3被视为测量数据，如图4所示.材料的实际参数与预测值如表2所示.与材料的实际参数相比，识别参数是非常准确的.

图4 正交异向简支梁的节点排列Fig.4 Nodal arrangements for the orthotropic simply beam

实例参数/(10-10Pa-1)预测值实际值预测值/实际值确定值迭代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] 孙秀山,黄立新,刘应华,等. 二维正交各向异性结构弹塑性问题的边界元分析[J]. 复合材料学报, 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] 黄立新，向志海，孙秀山，等.正交各向异性复合材料孔板性能参数识别测点的最优布置[J]. 工程力学, 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] 袁亚湘,孙文瑜. 最优化理论与方法[M]. 北京：科学出版社,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

湖南省自然科学基金资助项目(2017JJ2065)；教育部人文社会科学研究青年基金资助项目(10YJC630338)

TU313.3

A

1000-2537(2017)04-0062-06

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