第１２卷第４期 ２Ｏ１１年８月 理工大学学报（自然科学版） Ｊｏｕｒｎａｌ ｏｆ ＰＬＡ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｓｃｉｅｎｃｅ ａｎｄ Ｔｅｃｈｎｏｌｏｇｙ（Ｎａｔｕｒａｌ Ｓｃｉｅｎｃｅ Ｅｄｉｔｉｏｎ） Ｖｏ１．１２ Ｎｏ．４ Ａｕｇ．２０１ １ 基于线法的功能梯度材料断裂分析 燕秀发 。， 钱七虎 ， 方国强。， 赵敏福 ， 郭延宝 （１．理工大学工程兵工程学院，江苏南京２１０００７；２．９１５５０，辽宁大连ｌ１６０２３； ３．海军装备部，北京１００８４１；４．皖西学院，安徽六安２３７３００；５．装甲兵工程学院，北京１０００７２） 摘 要：为了克服材料非均匀性引起的数值困难，一种半解析数值方法——线法，被引入功能梯度材料的断 裂分析。通过有限差分将问题的控制方程半离散为定义在沿梯度方向离散节线上的常微分方程组，然后应 用Ｂ样条高斯配点法求解该常微分方程组。为了演示线法功能梯度材料断裂分析的方法，给出了指数梯度 含裂纹功能梯度材料板分别在恒定位移、弯曲载荷作用下应力强度因子的数值算例，同时给出了节线和配点 局部加密、区间映射以及变间距离散等提高计算精度和效率的线法断裂分析技巧。与相关问题理论解的对 比分析表明，该法的计算结果具有很高的精度。 关键词：功能梯度材料；线法；非均匀性；断裂分析；半解析法 中图分类号：Ｏ３４６．２ 文献标识码：Ａ 文章编号：１００９—３４４３（２０１０）０４—０３４６—０８ Ｍｅｔｈｏｄ ｏｆ ｌ ｉｎｅｓ ｆｏｒ ｆｒａｃｔｕｒｅ ａｎａｌｙｓｉｓ ｏｆ ｆｕｎｃｔｉｏｎａｌ ｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌｓ ＹＡＮ Ｘｉｕ—ｆａ ’ ， ＱＩＡＮ Ｑｉ—ｈｕ ， ＦＡＮＧ Ｇｕｏ—ｑｉａｎｇ。， ＺＨＡＯ Ｍｉｎ一－厂“ ， ＧＵＯ Ｙａｎ—ｂａｏ （１．Ｅｎｇｉｎｅｅｒｉｎｇ Ｉｎｓｔｉｔｕｔｅ ｏｆ Ｃｏｒｐｓ ｏｆ Ｅｎｇｉｎｅｅｒｓ，ＰＬＡ Ｕｎｉｖ．ｏｆ Ｓｅｉ．＆Ｔｅｅｈ．，Ｎａｎｊｉｎｇ ２１０００７，Ｃｈｉｎａ； ２．Ｕｎｉｔ Ｎｏ．９１５５０ ｏｆ ＰＬＡ，Ｄａｌｉａｎ １１６０２３，Ｃｈｉｎａ；３．Ｎａｖｙ Ｅｑｕｉｐｍｅｎｔ Ｄｅｐａｒｔｍｅｎｔ，Ｂｅｉｊｉｎｇ １００８４１，Ｃｈｉｎａ： ４．Ｗｅｓｔ Ａｎｈｕｉ Ｕｎｉｖｅｒｓｉｔｙ，Ｌｉｕａｎ ２３７３００，Ｃｈｉｎａ；５．ＰＬＡ Ａｒｍｏｒｉｎｇ Ｆｏｒｃｅ Ｅｎｇｉｎｅｅｒｉｎｇ Ｃｏｌｌｅｇｅ，Ｂｅｉｊｉｎｇ １０００７２，Ｃｈｉｎａ） Ａｂｓｔｒａｃｔ：Ｔｏ ｏｖｅｒｃｏｍｅ ｔｈｅ ｎｕｍｅｒｉｃａ１ ｄｉｆｆｉｃｕｌｔｉｅｓ ｃａｕｓｅｄ ｂｙ ｍａｔｅｒｉａｌ ｎｏｎｈｏｍｏｇｅｎｅｉｔｙ，ａ ｓｅｍｉ—ａｎａｌｙｔｉｃａｌ ｎｕ— ｍｅｒｉｃａ１ ｍｅｔｈｏｄ，ｍｅｔｈｏｄ ｏｆ １ｉｎｅｓ（ＭＯＬｓ）．ｗａｓ ｉｎｔｒｏｄｕｃｅｄ ｉｎｔｏ ｔｈｅ ｆｒａｃｔｕｒｅ ａｎａｌｙｓｉｓ ｏｆ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌｓ（ＦＧＭｓ）．Ｔｈｅ ｂａｓｉｃ ｉｄｅａ ｏｆ ｔｈａｔ ｍｅｔｈｏｄ ｗａｓ ｔｏ ｓｅｍｉ—ｄｉｓｃｒｅｔｉｚｅ ｔｈｅ ｇｏｖｅｒｎｉｎｇ ｅｑｕａｔｉｏｎｓ ｉｎｔｏ ａ ｓｙｓ— ｔｅｍ ｏｆ ｏｒｄｉｎａｒｙ ｄｉｆｆｅｒｅｎｔｉａｌ ｅｑｕａｔｉｏｎｓ（ＯＤＥｓ）ｂｙ ｍｅａｎｓ ｏｆ ｆｉｎｉｔｅ ｄｉｆｆｅｒｅｎｃｅ ａｐｐｒｏａｃｈ．Ｔｈｅ ＯＤＥｓ ｗｅｒｅ ｄｅ— ｆｉｎｅｄ ｏｎ ｔｈｅ ｄｉｓｃｒｅｔｅ ｌｉｎｅｓ ａｌｏｎｇ ｔｈｅ ｇｒａｄａｔｉｏｎ ｄｉｒｅｃｔｉｏｎ．Ｂｙ ｕｓｉｎｇ Ｂ ｓｐｌｉｎｅ ｃｏｌｌｏｃａｔｉｏｎ ｍｅｔｈｏｄ ａｔ Ｇａｕｓｓｉａｎ ｐｏｉｎｔｓ，ｓｏｌｕｔｉｏｎｓ ｔｏ ｔｈｅ ｐｒｏｂｌｅｍ ｗｅｒｅ ｏｂｔａｉｎｅｄ．Ｔｈｅ ｃｏｍｐｕｔａｔｉｏｎａｌ ｅｘａｍｐｌｅｓ ｏｆ ｔｈｅ ｓｔｒｅｓｓ ｉｎｔｅｎｓｉｔｙ ｆａｃｔｏｒｓ ｆｏｒ ｔｈｅ ｃｒａｃｋｅｄ ｐｌａｔｅ ｏｆ ＦＧＭｓ ｗｉｔｈ ａｎ ｅｘｐｏｎｅｎｔｉａｌ ｇｒａｄａｔｉｏｎ ｕｎｄｅｒ ｃｏｎｓｔａｎｔ ｄｉｓｐｌａｃｅｍｅｎｔ ａｎｄ ｂｅｎｄｉｎｇ ｌｏａｄ— ｉｎｇ ｗｅｒｅ ｐｒｅｓｅｎｔｅｄ ｔｏ ｄｅｍｏｎｓｔｒａｔｅ ｔｈｅ ｆｒａｃｔｕｒｅ ａｎａｌｙｓｉｓ ｍｅｔｈｏｄ ｏｆ ＦＧＭｓ．Ｉｎ ｓｕｃｈ ｅｘａｍｐｌｅｓ，ｓｏｍｅ ｈｕｍｅｒｉ— ｃａｌ ｔｅｃｈｎｉｑｕｅｓ，ｉｎｃｌｕｄｉｎｇ ｌｏｃａｌ ｒｅｆｉｎｅｍｅｎｔ ｏｆ ｄｉｓｃｒｅｔｅ ｌｉｎｅｓ ａｎｄ ｃｏｌｌｏｃａｔｉｏｎ ｐｏｉｎｔｓ，ｉｎｔｅｒｖａｌ ｍａｐ ａｎｄ ｖａｒｙｉｎｇ ｓｐａｃｅ ｄｉｓｃｒｅｔｉｚａｔｉｏｎ，ｗｅｒｅ ｅｍｐｌｏｙｅｄ ｔｏ ｉｍｐｒｏｖｅ ｔｈｅ ｃｏｍｐｕｔａｔｉｏｎａｌ ａｃｃｕｒａｃｙ ａｎｄ ｔｈｅ ｅｆｆｉｃｉｅｎｃｙ．Ｎｕｍｅｒｉｃａ１ ｒｅｓｕｌｔｓ ｏｆ ｔｈｅ ｅｘａｍｐｌｅｓ ｂｙ ＭＯＬｓ ａｒｅ ｑｕｉｔｅ ａｃｃｕｒａｔｅ ｃｏｍｐａｒｅｄ ｗｉｔｈ ｔｈｅ ｔｈｅｏｒｅｔｉｃａｌ ｓｏｌｕｔｉｏｎｓ ｔｏ ｔｈｅ ｃｏｒｒｅ— ｓｐｏｎｄｉｎｇ ｐｒｏｂｌｅｍｓ． Ｋｅｙ ｗｏｒｄｓ：ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌｓ；ｍｅｔｈｏｄ ｏｆ ｌｉｎｅｓ；ｎｏｎｈｏｍｏｇｅｎｅｉｔｙ；ｆａｃｔｕｒｅ ａｎａｌｙｓｉｓ；ｓｅｍｉ—ａｎａ— ｌｙｔｉｃａｌ ｍｅｔｈｏｄ 收稿日期：２００９—０１－０８． 基金项目：国家自然科学创新研究群体科学基金资助项目 （５１０２ｉ００１）；中国博士后科学基金资助项目 （２００８０４３１３４４）． 功能梯度材料是新一代非均匀材料，具有随空 间坐标变化的材料成分和特性，其独特之处在于可 作者简介：燕秀发（１９７２一），男，工程师；研究方向：防护工程数 值模拟与仿真；Ｅ—ｍａｉｌ：ｙｘｆｒｐ＠ｙａｈｏｏ．ＣＯｌＩ１．ｃｎ． 以通过定制材料成分和梯度来满足最终的需要，因 此在航空航天（热障涂层）、军事（梯度装甲）、光学 第４期 燕秀发，等：基于线法的功能梯度材料断裂分析 ３４７ （光纤）、核能（核反应堆第一壁）等领域都得到了广 泛的应用［１］。由于制备工艺的局限性和恶劣的工作 环境，断裂破坏是功能梯度材料的主要失效模式 ］。 因为功能梯度材料具有空间变化的材料属性参数， 使得描述该材料断裂问题的控制方程都成为变系数 偏微分方程，获得问题的解析解非常困难，所以数值 方法对于功能梯度材料的断裂分析具有特别重要的 意义。 目前，在功能梯度材料断裂分析中得到广泛应用 的是有限单元法和无网格法。根据实现材料物性函 数模拟的不同途径，功能梯度材料断裂分析的有限单 元法可分为２种类型，即均匀元法和等参梯度元法 （多重等参元法）。应用均匀元法，张幸红等［。 研究了 功能梯度材料裂纹应力强度因子与材料非均匀性参 数的关系。Ｇ．Ａｎａｌｓ等［４］分析了裂纹尖端应力场的 适用范围。Ｌｉ Ｃｈｕｎ－ｙｕ等［５］应用等参梯度元法计算 了该类材料环形裂纹的应力强度因子。Ｊ．Ｈ．Ｋｉｍ 等［６ 进一步将该法推广应用于正交各向异性功能梯 度材料结构的断裂分析。文献［７，８］应用无网格法计 算了功能梯度材料板边缘裂纹的应力强度因子，与解 析解的对比表明该计算结果具有较高的精度。 然而，在随后的研究中发现，上述方法存在一些 问题。均匀元法的基本原理是在单元内部采用零次 材料属性插值，通过不同单元问材料属性参数的变 化以分段近似来模拟功能梯度材料的材料梯度，其 根本缺陷是不能反映单元内部材料属性的变化，计 算精度依赖于网格细化程度［ｇ］。等参梯度元法扩展 了均匀材料有限单元法中等参变换的概念，应用材 料属性参数与位移场相同的插值函数进行模拟，因 此能够反映单元内部材料属性的变化，从而改善计 算精度。但是由于材料属性参数的连续变化可能导 致单元特性恶化，等参梯度元法的计算结果具有单 元类型敏感性和不确定性，表现为同一问题应用不 同的单元类型可能得到差异很大的计算结果ｒｌ 。 无网格法的问题在于，节点的影响域（支撑域）是重 叠的，与有限单元法中求解域被离散成有限个相互 连接的单元相比，无网格法的计算量远大于有限元 法的计算量，同时对材料属性参数的模拟进一步提 高了计算强度。特别是对于涂层、裂纹等属性参数 或待求场量具有高梯度性质的问题，计算强度的增 大更加明显，不能适应分析各种梯度形式对结构物 理特性影响的计算需要。 上述方法在功能梯度材料断裂分析中出现困难 的主要原因是，这些方法都是基于结构完全离散的 数值方法，计算误差不仅来源于求解区域的离散，也 来源于沿梯度方向材料属性参数的离散，连续变 化的材料属性参数导致了刚度矩阵数值积分的复 杂化以及可能的单元特性恶化。因此，李永等Ⅲ 提出了求解功能梯度材料板弯曲问题的半解析 法——康托洛维奇宏细观精化法（新康法）。与实 验结果的对比表明，该法具有很高的计算精度，且 便于分析各种参数对结构物理特性的影响，提供 了发展功能梯度材料结构分析半解析数值方法的 新思路。然而新康法针对特定问题选取特别的待 定函数形式，还不适合推广应用于一般性的功能 梯度材料结构分析问题。本文尝试将一种典型的 具有普遍意义的半解析数值方法——线法（ｍｅｔｈ— ｏｄ ｏｆ ｌｉｎｅｓ）引入功能梯度材料的断裂分析，通过具 体的算例演示该法的数值过程和计算技巧，并通 过与解析解的对比分析，最终表明线法在功能梯 度材料断裂分析中的有效性和适用性。 １ 问题描述 应力强度因子是断裂力学的重要参数，该参数 的计算是功能梯度材料断裂分析的重要内容。能否 克服裂纹尖端应力奇异性及其附近应力和应变的高 梯度，方便有效地得到应力强度因子的精确结果，是 检验断裂分析方法性能的有效手段。因此，本文要 研究的内容是沿Ｙ向梯度含裂纹功能梯度材料板应 力强度因子的求解问题，如图１所示。为使问题具 有代表性，分别考虑位移（图ｌ（ａ））和应力（图１（ｂ）） ２种载荷情况。 １０ ＲｌＩ ２／． ＋（ａ）受拉伸荷载作用 工 ／＋ Ｉ 。２ Ｒ／． ３＋ ｌ 。 （ｂ）受弯曲荷载作用 图１含裂纹功能梯度材料板 Ｆｉｇ．１ Ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌ ｐｌａｔｅ ｗｉｔｈ ａ ｃｒａｃｋ ３４８ 理工大学学报（自然科学版） 第１２卷 （１）均匀拉伸 ￡＝＝＝“ｆ０’ 为作用在板两侧边界处的拉伸位移。 （２）纯弯曲 一ａ￣［－Ｉ一 （ ＋ｓ）］。 式中： 为作用在板两侧边界处的弯曲应力，ｆｆｂｏ为 弯曲应力的幅值，ｗ为板的宽度，Ｓ为裂纹长度，Ｒ 为裂纹尖至板上端面的距离。 为与解析解比较，设材料弹性属性形式为 Ｅ—Ｅ（ ）， （１） ／ａ—Ｃ。 （２） 式中：Ｅ为弹性模量， 为泊松比，Ｃ为常数。考虑问 题的对称性，则相应载荷类型下待求问题的等效计 算构形如图２所示。 ｌ ｌ “ｒ０ ７ ７ Ｄ Ｄ ＼ Ｌ ￡ （ａ）受拉伸荷载作用（ｂ）受弯曲荷载作用 图２含裂纹功能梯度材料板的等效构形 Ｆｉｇ．２ Ｅｑｕｉｖａｌｅｎｔ ｍｏｄｅｌ ｏｆ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌ ｐｌａｔｅ ｗｉｔｈ ａ ｃｒａｃｋ ２ 数值实现 ２．１ 线法功能梯度材料断裂分析的控制方程 类似其他结构分析方法，线法也以位移作为基 本求解量，因此必须导出功能梯度材料用位移表达 的平衡方程和边界条件。 ２．１．１ 平衡方程 根据功能梯度材料制备完成后具有连续变化的 微结构和材料属性的特点，可将其处理为物性参数 是空间坐标连续函数的非均匀材料，而且功能梯度 材料的属性一般从一点到另一点逐渐变化，可以假 定给定一点的材料性质在任意方向上都是相同的。 因此，在连续介质意义上功能梯度材料被认为是各 向同性的非均匀实体。 综上所述，考虑功能梯度材料为物性参数是空 间坐标连续函数的各向同性非均匀材料，记 一 一（ｚ Ｙ ｚ），ｌＩ—ｌＩ ：（ Ｗ），ｉ一１，２，３。 为不失一般性，设材料参数为Ｅ＝Ｅ（ ）， ＝＝＝ Ｆ （ ），Ｇ—Ｇ（ ）一百７ ，Ｈ—Ｈ（ ）一 万 弋，则由物理方程、几何方程和应力平衡微分 １一 “ 方程可得功能梯度材料位移表示的平衡方程 Ｇ（ ＋ｕｓ，ｉｊ）＋２ＧＨ ＋（ ，＋筹 ）＋ （２ Ｈ＋２Ｇ ） 一０， （３） ｄ．２７ｉ ｄＸｉ 式中，ｉ、 ＝１，２，３。 ２．１．２边界条件 由物理方程、几何方程可得在应力边界Ｓ 上， Ｔｆ：［Ｇ（ ＋“ ）＋２ＧＨｕ ｋ， Ｉｎ』。 （４） 式中： 为作用在ｓ 上的表面力；ｎ』＝ＣＯＳ（Ｙ／， ）， 为边界面的外法线方向； 为Ｋｒｏｎｅｃｋｅｒ符号。 在位移边界Ｓ 上有 “ 一 ， （５） 式中， 为作用在Ｓ 上的位移。 由式（３）～（５），对比均匀材料位移形式的平衡 方程和边界条件可以看出，除位移边界条件外，功能 材料的控制方程均为包含材料属性函数的变系数偏 微分方程。 ２．２ 求解区域的离散和定义在节线上的控制方程 线法功能梯度材料结构分析的求解目标是定 义在沿梯度方向离散节线上的常微分边值问题， 首先必须将求解区域用若干沿梯度方向的节线离 散。考虑求解问题的数值特点，如图３所示，原问 题的等效计算构形为疏密相间的１～ｎ 条节线所 离散。图中ｉ为节线编号，并记各疏密区域分界节 线编号为 、 。、 、 ，求解区域边界节线编号为 １、 、 、嘶。ｈｉ为节线间距，前后相邻节线的间距 分别为ｈ斗 和ｈ卜 ，并记Ｈ。和Ｈ。分别为等距离散 区域（１）（４）和区域（３）（６）的节线间距。图中边界 节线外部编号为０、 、 ：、ｎ 的节线是引入边界条 件所需的虚拟线。 ０ ！’ Ｊ ＾ Ｉｒ． １ ｈ ｌ ｈ ｈ 十Ｊ ＪＨ２ Ｈ １ ２ ３ １Ｉ‘ ｎ】 ｆ一】 ｆ＋ｌ ｎ２ ３ ｌ ｎ ｒ 了 ’！ ｒ ｒ １ １ ｎ Ｉ ．．． ｎ５ ６ ｉ １ Ｔ （４） 区 （５） 区 ６ 图３等效计算构形中求解区域的离散 Ｆｉｇ．３ Ｄｏｍａｉｎ ｄｉｓｃｒｅｔｉｚａｔｉｏｎ ｏｆ ｔｈｅ ｐｒｏｂｌｅｍ ｉｎ ｔｈｅ ｅｑｕｉｖａｌｅｎｔ ｍｏｄｅｌ 第４期 燕秀发，等：基于线法的功能梯度材料断裂分析 ｃ ＋ 一ｃ ＋ Ｏ＇Ｏ ，３４９ 断裂问题是典型的高梯度问题，裂纹尖端附近 应力、应变具有剧烈变化的特点，同时该处还存在材 ｃ 。 料属性参数的变化。为兼顾计算精度和效率，求解 区域采取局部节线细化，节线细化区与粗糙区变间 （ ＋ 一（ Ｏ ｒ＋ ａｕ）Ｉ（ｘｊ，０） ’（１１）３，距过渡的离散方法，具体说明如下： （１）区别于裂纹边界和非裂纹边界，以裂纹尖 端为界，应用上下两排离散线； （２）在裂纹尖端附近，如图３中的（１）（４）区，设 式中，Ｊ— 。＋ｉ。 （４）在节线ｉ∈［ ， ］的端点（ ，一Ｓ）处有 （筹＋ ｄａｚｕ）（ａｅｉ，－Ｓ）一。， （１２） 置高度细化的等间距节线，以克服该处待求场量的 高梯度困难； （３）设置变间距（过渡）节线区，如图３中的（２） （５）区，但注意避免相邻节线间距比值过大； （４）在边界附近区域，如图３中的（３）（６）区，应 用相对粗糙的等间距节线，同样注意避免该区域离 散节线间距相对于过渡区邻近节线离散间距的较大 波动。 ２．２．１定义在离散节线上的控制方程 由式（１）～（３）可得图１、２中待求问题平面应力 情况下的控制方程，则相应的定义在节线ｉ上的平 衡方程为： ａ２Ｕ ｏ￣２Ｕ ．ｎ ．ｂ（ ｄｌ［（ ＋ （ ） ］＝０， （６） 争 ＋ ｄ。（ ） 一０。 （７） 式中：。一 ，６：＝＝ 一 一罟， 一 。其中， ＝ 。 ２．２．２ 离散节线端点处的边界条件 记ｚ 为节线 的 轴坐标，下标ｉ为节线编号。 由式（４）（５）有： （１）在节线ｉ一１的端点（０，Ｒ）处由约束条件 （避免结构发生平动）有： Ｉ（ｏ。Ｒ）一０， Ｉ（ｏ，Ｒ）一０。 （２）在节线ｉ∈［２， 。］的端点（ｚ ，Ｒ）处有： （Ｏｒｄｖ＋ 笔） Ｒ）＝＝：ｏ， （８） （Ｏｕｄ＋ ） 一０。 （９） ｖ（３）在节线ｉ∈［１， ］的端点（ ，Ｏ）处由位移 和应力连续性条件有： Ｉ（ ｆ，０）＝＝＝“Ｉ（ ，ｏ）’ Ｉ（ ，ｏ）＝＝： ｌ（ｚ，，０）， Ｏ ｕ ４Ｏｒ、ａａ）一，。。 （１３） ｚ（ｘｌ，－－Ｓ）２．３常微分边值问题的导出 ２．３．１ 区间映射 由图３中求解区域的离散方法和§２．２．２中离 散节线的端点边界条件可见，沿梯度Ｙ方向，上下两 组离散节线分别定义在不同的区间，即线１～ｎ。沿 Ｙ方向的定义区间为Ｅｏ，Ｒ］，线 ～ 沿该方向的 定义区间为［～Ｓ，Ｏ］。这种不规则的节线定义区 间，导致问题成为多点边值问题，非常不利于程序的 实现和求解效率的提高。因此需要应用区间映射技 术，将定义在不同区间的节线控制方程变换到同一 区间，使多点边值问题转变为两点边值问题。采用 如下的坐标变换： ｅ～Ｙ＿’ １ｌ （１４） 百一一百 Ｒ ≤Ｒ瓜口。Ｊ 葚 ，一 ≤－Ｓ一。 ≤。。ｏ｝Ｊ （１５） 则上、下两排节线沿梯度方向均被定义于［Ｏ，１］区 间，各离散节线端点坐标也成为（ｚ ，０）和（ｚｉ，１）。 ２．３．２ 节线位移函数对离散坐标偏导数的有限差 分近似 为将定义在节线上的控制方程转变为常微分方 程，应用由Ｔａｙｌｏｒ公式得到的具有二阶精度的非等 距三点中心差分近似该方程中待求场量对离散坐标 ｚ的偏导数，记节线位移函数 一 （ ）＝ｕ（ｘ ， ）， 一 （ ）一 （ｚ ， ）及塑ｄ７３ｉｄｙ一 ， 一 ， 鲁 ： ；，，孥一 ，则对于图３中节线间距为 ｈ斗 、ｈｉ、ｈＨ的区域（２）（５）有： ３５Ｏ 理工大学学报（自然科学版） 第ｌ２卷 （嘉）。 ２Ｅｈ ｌ一（＾ ＋ｈ １） ｆ＋ｈ １ 斗１］ ｈｉ １ｈ，（ｈｌ－１＋ｈ ） （舅） ｒａｚ￡、 一雎１ 斗１＋（＾ ～硅１）“ｉ一＾ ｚ‘卜１］ 几一——— —一’ ｒ３ｖ、几 一＾ １—— 斗１＋（— ～＾ｚ＿１） 一 —一卜１］ ’ ，３ 、一 ｌ ｚ￡ １＋（＾ 一 ｚＩ１） ：一＾ Ｍ：～１］ 一——— —一’ ｒ ａ 、 硅１ ｌ＋（＾ 一 １） ：一ｈｆ３７￣＇１］ 、ａｚａ ｈｉｈ卜１（＾卜１＋ｈ ） 。 （１６） 如果式（１６）中ｈ斗 一ｈ 一ｈ 。一Ｈ 或Ｈ。，则 得到适用于图３中的（１）（３）（４）（６）等区域的等距差 分公式。当ｉ∈（１，哟）和ｉ∈（ ， ）时，将式（１６） 代入式（６）（７），则定义在节线上的控制方程由偏微 分方程组转变为以节线位移函数为未知量的常微分 方程组。显然，对于边界节线即图３中的线１、 。和线 、 ，为了应用式（１６），针对本文问题需要得到虚 拟线上的 。、 、 ：、 、 、 、 ：、口 、Ｕ ｏ、ｕ：ｏ、 、 等辅助量，这些辅助量通过引入该线处的边 界条件获得。 ２．３．３边界线上边界条件的引入和虚拟线上的辅 助量 （１）拉伸和弯曲２种载荷情况下虚拟线ｉ＝０ 上的辅助量７３ｏ、 。 、，．．．．． ．，。．．．． 、，．．．．．．．．． 根据问题的对称性，在边界线（ｉ一１）处有边 界条件： “１ Ｏ， ］ Ｇ［（ ）ｄ～ ＋（ ）ｄｚ ］一０。）｝ 由式（１６）的第４式有（ ） 一 。 代入式（１７）得： ３７ｏ ＝＝＝ ’Ｉ３７２。】 （１８） （２）拉伸和弯曲２种载荷情况下虚拟线 ： ： 上的辅助量 “ 、 。 在ｉ一 处有边界条件： ’ Ｇ［（ ＋（ ］一０。 ｊ 由式（１６）第３、４式有： 一 ， ｃ 一 。 将式（２Ｏ）代入式（１９），并注意式（１）（２），可得： Ｕｎ４＋１十 Ｈ ，１ 一 ｎ４＋ｌ十 Ｈ ，ｌ ＞ “ 一“ ４＋１十 Ｈ１ＩＺ－， ￣ ，ｌ （２１） 口 一３７ｎ４＋１十 Ｈ。 。Ｊ 类似 、 的求解过程可得到均匀拉伸载荷下 虚拟线ｉ＝ ：和ｉ： ；上的辅助量： ， （２２） ｎ！ ， ”！ 类似 ：、 、 ：的求解过程可得到弯曲载 荷下虚拟线ｉ：ｎ 和ｉ：ｎ 上的辅助量： ｏ— 一２Ｈ。 ＋ ［１一 （ ＋ｓ）］’］ 一 一２Ｈ。 。， Ｊ “ 一“ 一 一２Ｈ ：。一下２Ｈ２ａ￣ 一生警 ＋｝ ＋ｓ）， Ｉ咄一 。一。一２Ｈ： 。。 ｊ ｒ９ 、 ０一 一２Ｈ 十 Ｅ１一 （ ＋ｓ）］，１ 一 一２Ｈ：乱 ， Ｉ“ 一 一 一２Ｈｚ 一下２Ｈ２ａ￣Ｅ 一４面Ｈ２ａｂｏ＋｝ ＋ｓ）， ｝ ＝ ： 一 一２Ｈ 。 ｊ （２４） 式中， Ｅ／（卜 一￡（ ）＝ 。 ２．３．４定义在节线上的常微分方程 应用式（１６）近似式（６）（７）中位移函数对离散坐 标ｚ的偏导数，对边界线注意式（１８）（２１）（２２）（拉伸 载荷）式（２３）（２４）（弯曲载荷），然后利用式（１４）进行 区间变换，则可得定义在节线ｉ∈［１， 。］上的原问 题控制方程的常微分方程近似式： ＝Ｒ ＋Ｒｚ豆 ＋Ｒ。 。 （２５） 式中：一Ａ，Ｂ一和 为系数矩阵； 一 ｄＵ，对于均匀拉 伸载荷待求位移函数 一（ ｌ Ｓｔ２ ２…“ … 一１ 第４期 燕秀发，等：基于线法的功能梯度材料断裂分析 ７Ｊ 一。 ） ，对于弯曲载荷待求位移函数 ：（ “： ７）２… …Ｕｎ３－１ ７）ｎ３－ｌ ｎ３ Ｈ３）Ｔ。 利用式（１５），经过与式（２５）类似的过程可得定 义在节线ｉ∈［ ， ］上的原问题控制方程的常微 分方程近似式： 一一ＳＡ Ｕ ＋Ｓ Ｂ Ｕ＋Ｓ２Ｆ， （２６） 式中：～Ａ、里和Ｅ为系数矩阵； ＝ ，对于均匀拉 ｕｑ－ 伸载荷待求位移函数， 一（ “ ＋－＂Ｏｎ ＋－…甜， … ，～。 ，一。 ） ，对于弯曲载荷待求位移函数 ＝（ ４ ４ ～＋１ ４＋１…ｚ￡ 』… ７ ７） 。 ２．３．５端点边界条件的常微分形式 将式（１４）代入式（８）～（１３），并注意对于ｉ∈ ［２，嘞一１３和ｉ∈［ ＋１， 一１３等内部节线，应用 式（１６）中的第３、４式近似式（８）～（１３）中的位移函 数对离散坐标 的偏导数。但对于边界节线，为获 得与内部节线相同的近似精度，需应用由Ｎｅｗｔｏｎ 等距插值公式得到的式（２７）来近似各边界节线端点 边界条件方程中对位移坐标ｚ的偏导数ｍ］。 （ ） ： ，１ｌ （２７） ，ａ厂、一３ａ 一 一４ｆ￣２＾ －１＋ｆｍ－２ ｌ。Ｊ 式中：，为待求位移函数 或 ．ｆ、ｍ分别为求解域 两侧边界节线的编号，ｚ一１、 ４，ｍ一现、 ；ｈｆ— Ｈ ，ｈ 一Ｈ：。于是可得ｉ∈ＥＬ， 。］各节线端点边界 条件的常微分方程近似式 Ｃ Ｉ ｕ ＋ｔ５０ Ｊ ｒ，Ｕ＋ＲＧ Ｉ ｒｋ一０。 （２８） 式中：ｋ一０，１，Ｆ０一（ ０），Ｆ１＝（ｚ １）。 利用式（１５），经过与式（２８）类似的过程可得 ｉ∈［ ， ］各节线端点边界条件的常微分方程近 似式 Ｉ【， 一Ｓ Ｉ 一Ｓ Ｉ ｒｋ一０。 （２９） 经过上述区间映射和有限差分过程，原求解域内 的偏微分边值问题转变为定义在各离散节线上的由 式（２５）（２６）（２８）（２９）构成的常微分两点边值问题。 ３常微分边值问题的求解和应力强度 因子的计算 ３．１ 常微分边值问题的求解 将原求解域的偏微分边值问题转化为定义在离 散节线上的常微分边值问题后，常微分边值问题的 求解成为关键。在求解常微分边值问题的方法中， 打靶法和样条函数配点法是较为优越的算法。然 而，Ｊｏｎｅｓ等的研究表明，应用打靶法求解线法常微 分方程组时存在计算结果不稳定的可能［１３３。基于 这一考虑，本文选择样条函数配点法中的Ｂ样条 Ｇａｕｓｓ配点法求解线法常微分边值问题。 Ｂ样条Ｇａｕｓｓ配点法求解常微分边值问题的基 本原理是：以Ｂ样条基样条插值函数作为问题的近似 解，通过使其在各配点区间的Ｇａｕｓｓ映射点和边界点 处满足常微分方程和边界条件确定基函数的系数，由 此得到近似解。求解节线常微分方程组得到离散节 线位移近似函数后，对位移近似函数求导则可得沿节 线方向的应变近似函数，插值可得节线上任意点在该 方向上的应变值，由相邻节线对应点位移数值的有限 差分，则可得离散方向上的应变数值。将所得应变值 代人物理方程，于是进一步得到该点处的应力数值。 ３．２计算应力强度因子的方法 裂纹尖端附近的应力场和位移场取决于应力强 度因子Ｋ，对功能梯度材料裂纹尖端场的研究表 明，该类材料裂纹尖端处应力场与均匀材料具有相 同的函数形式，但应力强度因子由裂纹尖端处的材 料参数决定，图１、２中所示裂纹应力强度因子的定 义 为： Ｋ 一 ０譬 。Ｖ ｒ 式中： 为裂纹面的张开位移，ｒ为该张开位移点到 裂纹尖端的距离， 为裂纹尖端处的弹性模量。将 用 ／（１一／ｚ。）代替可得平面应变情况下应力强 度因子的定义。 基于位移的外推法是断裂问题中计算应力强度 因子的基本方法。但是，尽管应用线法可以得到节 线上任一点的位移值，但在裂纹尖端ｒ＝＝＝０处，有限 差分不能真正反映该处位移对坐标导数趋于无穷的 性质。因此，在靠近裂纹尖端处线法计算结果的误 差较大，不能直接应用上述极限过程。基于上述考 虑，通过下面的方法计算应力强度因子，定义Ｋ 为 名义应力强度因子，则 Ｋ 一 。 以Ｋ 为纵坐标，ｒ为横坐标，将有关点绘出，用最小 二乘法处理，通过这些点绘出一条最佳拟合直线后， 该直线至纵坐标轴的交点即为是 的估算值。 ３．３进一步提高计算精度的方法 令图１、２中Ｖ ＝１，Ｌ＝：：３Ｗ，Ｒ＝０．８，Ｓ一 理工大学学报（自然科学版） 第１２卷 ０．２，“∞一１， 一１，并记Ｅｏ：Ｅ（Ｏ），Ｅ１一Ｅ（一Ｓ）， Ｅ 一Ｅ（Ｒ）。为了与已知的解析解［１印对比，考虑式 （１）（２）中的材料梯度模型为 Ｅ—Ｅｏｅ ． ｔｚ ０．３， 分别计算各载荷在Ｅ２／Ｅ ＝０．５、２两种梯度水平下 的应力强度因子，口和Ｅ。可根据材料内部弹性模量 空间变化情况的数据拟合得到。 考虑与已知解析解的前提一致，计算问题的平面 应变情况，只需令Ｅ—Ｅ／Ｃ１－－ｔｚ ）和／ｚ一 ／（１一 ）， 即可方便地将平面应力问题转换为平面应变问题。 计算过程中取协调一致的封闭单位体系，各计算参数 的量纲为相应的单位力、单位长度、单位弹性模量。 考虑到裂纹尖端附近存在着高梯度的应力和应 变并兼有材料属性梯度，为了获得良好的计算精度， 不仅需要裂纹尖端附近离散节线高度细化，还可以 根据Ｂ样条Ｇａｕｓｓ配点法求解常微分方程的数值 原理，采取进一步提高计算精度的措施。对于二阶 常微分方程，Ｂ样条Ｇａｕｓｓ配点法的精度可达到 ２ １曼 ０ ３ …Ｏ（Ａ ５ Ｚ＋ ｋ），２ 其中，△ ５ ｌ ５ 为配点区间的最大长度，志为各 Ｏ 配点区间的Ｇａｕｓｓ映射点数。因此，可以通过在求解 区间内设置较多的配点区间以减小△ （类似于有 限单元法中的Ｈ算法）和较多的Ｇａｕｓｓ映射点以增 大是（类似于有限单元法中的Ｐ算法），并特别注意配 点区间向裂纹尖端局部加密以最大限度降低该处的 计算误差，进一步提高断裂分析的计算精度。 ４应力强度因子的数值结果 在裂纹尖端附近位移偏导数趋向于无穷，且具 有高梯度变化的特点，所以在该处应用有限差分近 似位移偏导数存在较大误差，相应的数值结果在裂 纹尖端附近并不准确，这种情况同其他结构分析方 法如非奇异元有限单元法都是类似的。因此，应力 强度因子应由距裂纹尖端一定距离处的名义应力强 度因子通过最小二乘法得到的最佳拟合直线外推获 得。图４是均匀拉伸载荷情况下用§３．２中的方 法，由最佳拟合直线获得应力强度因子的外推过程。 图中拟合点为各配点区间的Ｇａｕｓｓ映射点，为与解 析解对比，各拟合点纵坐标为归一化后的名义应力 强度因子Ｋ ／（ 。 ７ｃｓ），其中， ：Ｅ “ ／ＥＬ（１一 。）］。外推直线与纵轴的交点即为归一化后的应力 强度因子。用类似的方法可得弯曲载荷情况下的归 一化应力强度因子Ｋ ／（ 加ｖｒ￣Ｓ）。 Ｉ ｏ Ｅ／Ｅ 卸．５ Ｉ 。 Ｉ＋ ／Ｅ，＝２．０ ｌ ０ ０．０５ ０．１０ ０．１５ ０．２０ ， 图４拉伸载荷下归一化的应力强度因子 Ｆｉｇ．４ Ｎｏｒｍａｌｉｚｅｄ ｓｔｒｅｓｓ ｉｎｔｅｎｓｉｔｙ｛ａｃｔｏｒｓ ｕｎｄｅｒ ｔｈｅ ｔｅｎｓｉｏｎ ｌｏａｄｉｎｇ 表１是图１、２中各载荷在两种梯度情况下应力 强度因子线法数值解与解析解［１ｓ］的对比。比较表 明，线法数值解与解析解符合得较好。同时表中结 果也说明，材料梯度变化情况直接影响计算精度。 因此，上述算例不仅证实了线法功能梯度材料断裂 分析的有效性，也进一步表明该法对于材料梯度参 数波动具有良好的适应性。分析线法的数值原理可 以发现：在非梯度方向，有限差分的离散精度是 ０（＾。）（＾为节线间距）；在梯度方向，节线常微分方 程的求解精度是０（△…２＋ｋ）。因此，加密节线和配点区 间或增加Ｇａｕｓｓ映射点数可以进一步提高精度，缩 小因材料梯度增大引起的计算误差。 表１ 应力强度因子对比 Ｔａｂ．１ Ｃｏｎｔｒａｓｔ ｏｆ ｓｔｒｅｓｓ ｉｎｔｅｎｓｉｔｙ ｆａｃｔｏｒｓ ５ 结 论 本文首先导出了功能梯度材料位移形式的平衡 方程和边界条件，然后以拉伸、弯曲等载荷作用下指 数梯度材料平面裂纹问题应力强度因子的计算为 例，演示了线法功能梯度材料断裂分析的基本过程。 该法的基本思想是沿着与材料梯度无关的方向将求 解区域半离散化，通过有限差分将问题的控制方程 转变成为定义在离散节线上的常微分方程组，应用 Ｂ样条函数Ｇａｕｓｓ配点法求解该常微分方程组得到 问题的解答。此外，在算例中同时给出了节线变间 距离散、区间映射以及节线和配点区间局部加密等 第４期 燕秀发，等：基于线法的功能梯度材料断裂分析 ３５３ 提高计算精度与效率的数值技巧。线法计算结果与 解析解的对比表明了线法在功能梯度材料断裂分析 中的有效性。 由线法的求解过程可见，与有限单元法和无网 格法等基于结构完全离散并求解节点处代数方程的 ［５］ＬＩ Ｃｈｕｎ－ｙｕ，ＺＯＵ Ｚｈｅｎｇ－ｚｈｕ．Ｓｔｒｅｓｓ ｉｎｔｅｎｓｉｔｙ ｆａｃｔｏｒ ｆｏｒ ａ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌ ｃｙｌｉｎｄｅｒ ｗｉｔｈ ａｎ ｅｘｔｅｒｎａｌ ｃｉｒ— ｃｕｎｄｅｒｅｎｔｉａｌ ｃｒａｃｋ［Ｊ］．Ｆａｔｉｇｕｅ＆Ｆｒａｃｔｕｒｅ ｏｆ Ｅｎｇｉｎｅｅｒｉｎｇ Ｍａｔｅｒｉａｌｓ ８Ｌ Ｓｃｉｅｎｃｅ，１９９８，２１（９）：１４４７—１４５７． ［６］ＫＩＭ Ｊ Ｈ，ＰＡＵＬＩＮＯ Ｇ Ｈ．Ｍｉｘｅｄ－ｍｏｄｅ ｆｒａｃｔｕｒｅ ｏｆ ｏｒ－ ｔｈｏｔｒｏｐｉｃ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌｓ ｕｓｉｎｇ ｇｒａｄｅｄ ｆｉｎｉｔｅ 完全数值方法不同，线法直接求解的是结构半离散 后，定义在沿梯度方向节线上保留了材料属性参数 变化信息的常微分方程，具有典型的半解析方法的 性质。因此，该法能够有效地反映功能梯度材料连 ｅｌｅｍｅｎｔｓ ａｎｄ ｔｈｅ ｍｏｄｉｆｉｅｄ ｃｒａｃｋ ｃｌｏｓｕｒｅ ｍｅｔｈｏｄ［Ｊ］．Ｅｎｇｉ— ｎｅｅｒｉｎｇ Ｆｒａｃｔｕｒｅ Ｍｅｃｈａｎｉｃｓ，２００２，６９（１４－１６）：１５５７－１５８６． ［７］陈建，吴林志，杜善义．采用无单元法计算含边沿裂 续变化的材料属性，不需要像均匀元法分段近似材 料属性函数，也不需要类似等参梯度元法进行附加 的材料属性参数的插值近似。对比新康法，线法对 问题的求解直接从控制方程出发，不需要事先选取 满足边界条件的特定待求场函数，求解的问题更加 广泛，是一种具有普遍意义的半解析数值方法。 进一步分析线法的数值特点可以看出，线法实 质上是在１个或２个方向（三维问题）用有限差分求 近似解，而在材料梯度方向用Ｂ样条Ｇａｕｓｓ配点法 求高精度解，类似于有限元分析的网格局部加密和 采用高阶位移模式，因此非常符合目前具有单向梯 度的大多数功能梯度材料的特点。此外，应用Ｂ样 条Ｇａｕｓｓ配点法求解常微分边值问题，可以通过缩 小配点区间长度或增加配点数，改善局部的计算精 度，使线法具有很高的应用灵活性，特别适合求解裂 纹、应力集中和材料梯度陡变等高梯度问题。 参考文献： ［１］黄旭涛，严密．功能梯度材料：回顾与展望［Ｊ］．材料 科学与工程，１９９７，１５（４）：３５－３８． ＨＵＡＮＧ Ｘｕ－ｔａｏ，ＹＡＮ Ｍｉ．Ｒｅｖｉｅｗ ａｎｄ ｐｒｏｓｐｅｃｔｓ ｏｆ ｆｕｎｃｔｉｏｎａｌ ｇｒａｄｉｅｎｔ ｍａｔｅｒｉａｌｓ［Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ｍａｔｅｒｉａｌｓ Ｓｃｉ－ ｅｎｃｅ ａｎｄ Ｅｎｇｉｎｅｅｒｉｎｇ，１９９７，１５（４）：３５—３８．（ｉｎ Ｃｈｉｎｅｓｅ）． ［２］ＬＡＮＵＴＴＩ ＪＪ．Ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌｓ：ｐｒｏｐｅｒ— ｔｉｅｓ，ｐｏｔｅｎｔｉａｌ ａｎｄ ｄｅｓｉｇｎ ｇｕｉｄｅｌｉｎｅｓ［Ｊ］．Ｃｏｍｐｏｓｉｔｅｓ Ｅｎｇｉｎｅｅｒｉｎｇ，１９９４，２４（４）：８１—９４． ［３］张幸红，李亚群，韩杰才．ＴｉＣ—Ｎｉ系功能梯度材料的断 裂力学分析ＥＪ］．复合材料学报，２００１，ｌ１（４）：８７—９２． ＺＨＡＮＧ Ｘｉｎｇ—ｈｏｎｇ．ＬＩ Ｙａ－ｑｕｎ，ＨＡＮ Ｊｉｅ－ｃａｉ．Ｃｒａｃｋ ｐｒｏｂｌｅｍ ａｎａｌｙｓｉｓ ｏｆ ＴｉＣ－Ｎｉ ＦＧＭ ｕｓｉｎｇ ｆｉｎｉｔｅ ｅｌｅｍｅｎｔ ｍｅｔｈｏｄ［Ｊ］．Ａｃｔａ Ｍｅｃｈａｎｉｃａｌ Ｃｏｍｐｏｓｉｔｅ Ｓｉｎｉｃａ，２００ １， ｌ１（４）：８７－９２．（ｉｎ Ｃｈｉｎｅｓｅ）． ［４］ ＡＮＡＬＳ Ｇ，ＬＡＭＢＲＯＳ Ｊ，ＳＡＮＴＡＲＥ Ｍ Ｈ．Ｄｏｍｉ— ｎａｎｃｅ ｏｆ ａｓｙｍｐｔｏｔｉｃ ｃｒａｃｋ ｔｉｐ ｆｉｅｌｄｓ ｉｎ ｅｌａｓｔｉｃ ｆｕｎｃｔｉｏｎ— ａｌｌｙ ｍａｔｅｒｉａｌｓ［Ｊ］．Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｆｒａｃｔｕｒｅ， ２００３，１１５（２）：１９３—２０４． 纹功能梯度材料板的应力强度因子［Ｊ］．工程力学， ２０００，１７（５）：１４０－１４４． ＣＨＥＮ Ｊｉａｎ，ＷＵ Ｌｉｎ－ｚｈｉ，ＤＵ Ｓｈａｎ－ｙｉ．Ｅｖａｌｕａｔｉｎｇ ＳＩＦ ｏｆ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｐｌａｔｅ ｗｉｔｈ ｆｉｎ ｅｄｇｅ ｃｒａｃｋ ｂｙ ｅｌｅｍｅｎｔ—ｆｒｅｅ ｍｅｔｈｏｄ［Ｊ］．Ｅｎｇｉｎｅｅｒｉｎｇ Ｍｅｃｈａｎｉｃｓ， ２０００，１７（５）：１４０－１４４．（ｉｎ Ｃｈｉｎｅｓｅ）． ［８］ＲＡＯ Ｂ Ｎ，ＲＡＨＭＡＮ Ｓ．Ｍｅｓｈ－ｆｒｅｅ ａｎａｌｙｓｉｓ ｏｆ ｃｒａｃｋｓ ｉｎ ｉｓｏｔｒｏｐｉｃ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌｓ［Ｊ］．Ｅｎｇｉ— ｎｅｅｒｉｎｇ Ｆｒａｃｔｕｒｅ Ｍｅｃｈａｎｉｃｓ，２００３，７０（１）：卜２７． ［９］ ＡＮＡＬＳ Ｇ，ＳＡＮＴＲＥ Ｍ Ｈ，ＬＡＭＢＲＯＳ Ｊ．Ｎｕｍｅｒｉｃａｌ ｃａｌｃｕｌａｔｉｏｎ ｏｆ ｓｔｒｅｓｓ ｉｎｔｅｎｓｉｔｙ ｆａｃｔｏｒｓ ｉｎ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌ［Ｊ］．Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｆｒａｃｔｕｒｅ， ２０００，１０４（２）：１３１－１４３． ［１Ｏ］ＫＩＭ Ｊ Ｈ，ＰＡＵＬＩＮＯ Ｇ Ｈ．Ｉｓｏｐａｒａｍｅｔｒｉｃ ｇｒａｄｅｄ ｆｉ— ｎｉｔｅ ｅｌｅｍｅｎｔｓ ｆｏｒ ｎｏｎｈｏｍｏｇｅｎｏｕｓ ｉｓｏｔｒｏｐｉｃ ａｎｄ ｏｒｔｈｏ— ｔｒｏｐｉｃ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌｓ［Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ａｐｐｌｉｅｄ Ｍｅｃｈａｎｉｃｓ，２００２，６９（３）：５０２—５ １４． ［１１］李永，宋健，张志民．非均质梯度材料复合结构的 Ｋａｎｔｏｒｏｖｉｃｈ宏细观精化解法［Ｊ］．中国科学（Ｅ辑）， ２００３，３３（１）：２９—４１． ＬＩ Ｙｏｎｇ，ＳＯＮＧ Ｊｉａｎ，ＺＨＡＮＧ Ｚｈｉ－ｍｉｎ．Ｋａｎｔｏｒｏｖｉｃｈ ｍｅｔｈｏｄ ｏｆ ｍａｃｒｏ ｏｒ ｍｅｓｏ—－ｓｃｏｐｉｃ ｓｏｌｕｔｉｏｎ ｆｏｒ ｎｏｈｏｍｏ。‘ ｇｅｎｅｏｕｓ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌｓ ｃｏｍｐｏｓｉｔｅ［Ｊ］． Ｊｏｕｒｎａｌ ｏｆ Ｃｈｉｎａ Ｓｃｉｅｎｃｅ（Ｓｅｒｉｅｓ Ｅ），２００３，３３（１）：２９— ４１．（ｉｎ Ｃｈｉｎｅｓｅ）． ［１２］凌复华，殷学纲，何冶奇．常微分方程数值方法在力学 中应用［Ｍ］＿重庆：重庆大学出版社，１９９０． ［１３］ＪＯＮＥＳ Ｄ Ｊ，Ｓ０ＵＴＨ Ｊ Ｃ，ＫＬＵＮＫＥＲ Ｅ Ｂ．Ｏｎ ｔｈｅ ｎｕｍｅｒｉｃａｌ ｓｏｌｕｔｉｏｎ ｏｆ ｅｌｌｉｐｔｉｃ ｐａｒｔｉａｌ ｄｉｆｆｅｒｅｎｔｉａｌ ｅｑｕａ— ｔｉｏｎ ｂｙ ｔｈｅ ｍｅｔｈｏｄ ｏｆ ｌｉｎｅｓ［Ｊ］．Ｊｏｕｒｎａｌ ｏｆ Ｃｏｍｐｕｔａ— ｔｉｏｎａｌ Ｐｈｙｓｉｃｓ，１９７２，９（２）：４９６—５２７． ｒ１４］ＲＡＶＩＣＨＡＮＤＲＡＮ Ｋ Ｓ，ＢＡＲＳＯＵＭ Ｉ．Ｄｅｔｅｒｍｉｎａｔｉｏｎ ｏｆ ｓｔｒｅｓｓ ｉｎｔｅｎｓｉｔｙ ｆａｃｔｏｒ ｆｏｒ ｃｒａｃｋｓ ｉｎ ｆｉｎｉｔｅ－ｗｉｄｔｈ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｍａｔｅｒｉａｌ［Ｊ］．Ｉｎｔｅｒｎａｔｉｏｎａｌ Ｊｏｕｒｎａｌ ｏｆ Ｆｒａｃｔｕｒｅ，２００３，１２１（２）：１８３—２０３． ［１５］ＥＲＤＯＧＡＮ Ｆ，ｗＵ Ｂ Ｈ．Ｔｈｅ ｓｕｒｆａｃｅ ｃｒａｃｋ ｐｒｏｂｌｅｍ ｆｏｒ ａ ｐｌａｔｅ ｗｉｔｈ ｆｕｎｃｔｉｏｎａｌｌｙ ｇｒａｄｅｄ ｐｒｏｐｅｒｔｉｅｓ［Ｊ］． Ｊｏｕｒｎａｌ ｏｆ Ａｐｐｌｉｅｄ Ｍｅｃｈａｎｉｃｓ，１９９７，６４（３）：４４９—４５６． （责任编辑：汤雪峰）
