微分求积法计算双材料反平面V形切口应力奇性指数

PAGEPAGE10基金项目：安徽省自然科学基金资助项目(1808085ME147)作者简介：葛仁余（1969-)，男，博士，副教授，研究方向：断裂力学和计算固体力学，Email：gerenyu@张佳宸(1996-)，男，安徽工程大学，硕士研究生，研究方向：计算固体力学，Email：381610972@微分求积法计算双材料反平面V形切口应力奇性指数摘要：文章运用微分求积法计算双材料反平面V形切口端部应力奇异性指数。首先基于弹性力学理论，将切口端部应力奇异性指数的计算转化为常微分方程组的特征值问题，再由微分求积法理论将常微分方程组的特征值问题转化为标准型广义代数方程组特征值问题，采用QR法可一次性地计算出双材料反平面V形切口端部应力奇异性指数及其特征角函数。数值计算结果表明：微分求积法对于不同材料参数组合进行了计算，可一次性解出双材料反平面V形切口端部所有的应力奇性指数及其相应的特征角函数，计算值与已有文献计算结果吻合良好，证明了本文数值计算方法的有效性和计算精度。关键词：反平面V形切口；双材料，应力奇异性，微分求积法Thecalculationofstresssingularityordersofanti-planeV-notchinbi-materialsusingDifferentialquadraturemethodAbstract:Inthispaper,thedifferentialquadraturemethodisusedtocalculatethesingularityordersofthestressattheendoftheBi-materialanti-planeV-notch.Firstly,basedonthetheoryofelasticity,thecalculationofthestresssingularityordersatthenotchendistransformedintotheeigenvalueproblemofordinarydifferentialequations,Secondly,theeigenvalueproblemofordinarydifferentialequationsistransformedintotheeigenvalueproblemofstandardgeneralizedalgebraicequationsbythetheoryofdifferentialquadrature,Then,QRmethodcanbeusedtocalculatethestresssingularityindexanditscharacteristicanglefunctionattheendofBi-materialanti-planeV-notch.Thenumericalresultsshowthat:Whenthedifferentialquadraturemethodisusedtocalculatethecombinationofdifferentmaterialparameters,thenallthestresssingularityordersattheendoftheV-notchanditscorrespondingcharacteristicanglefunctioncanbesolvedatonetime,andthecalculatedresultsareingoodagreementwiththoseintheliterature,whichprovestheeffectivenessandaccuracyofthenumericalmethod.Keywords:Anti-planeV-notches;Bimaterial;Stresssingularity;Differentialquadraturemethod1、引言近几十年来，人们一直在考虑不同材料的粘结楔块端部应力奇异现象，Bogy（1971）[1]和Dempsey、Sinclair（1979、1981）[2,3]研究了不同材料的粘结弹性楔块的平面应力集中问题。Ma和Hour（1989）[4]分析了各向同性和各向异性的异质双材料反平面V形切口尖端应力奇异性。Ting（1986）[5]研究了各向同性和各向异性材料在载荷作用下双材料界面裂纹尖端的平面应力奇异性。Kargarnovin等(1997)[6]和Shahani(1999)等[7]通过复势函数法配合有限Mellin变换求解问题的控制微分方程，获得了不同的边界条件下各向同性和各向异性反平面切口尖端位移场和应力场的解析解，以及反平面切口端部应力奇异性指数的显式表达式。平学成等[8]提出了一个基于位移的分析尖劈端部奇性位移场和应力场反平面问题的非协调元特征法，葛仁余等[9]采用插值矩阵法研究了幂次硬化材料反平面V形切口尖端塑性应力奇异性问题；葛大丽[10]采用插值矩阵法研究了各向同性双材料反平面V形切口尖端应力奇异性问题。钱俊等[11]先用特征函数法导出了反平面复合材料切口的特征方程，后采用Müller法计算了解的序列；PageauSS等[12]采用有限元特征法研究了平面各向异性复合材料切口在反平面载荷作用下的奇异应力状况。微分求积法是Bellman和Casti于1971年提出的基本理论[13]，由于快速计算机及仿真技术的发展而引起了广大学者们的广泛关注。该方法具有公式简单、编程方便、计算量少和精度高等优点而广泛应用于工程结构的动力学特性分析研究中。在国内本文首次将微分求积法应用到工程断裂力学的数值计算研究领域，首先根据弹性力学和断裂力学基本理论，将异质双材料反平面V形切口端部应力奇异性指数的计算转化为关于应力奇异性指数的常微分方程组的特征值问题，再基于微分求积法理论，将常微分方程组的特征值问题转化为标准型广义代数方程组特征值问题，采用QR法可一次性地计算出异质双材料反平面V形切口端部应力奇异性指数及其特征角函数。2、双材料反平面V形切口图1所示为双材料反平面V形切口，在切口尖端处同时定义一个的直角坐标系和一个的极坐标系，和为两个楔形尖劈的张角，和为两种材料的剪切弹性模量，和为切口边界，为两种材料的交界。对于一个反平面V形切口构件，根据线弹性理论分析，Champion等[14]将反平面V形切口尖端附近区域2个尖劈域内的的位移场表达成关于径向的一系列级数渐近展开形式：图1双材料反平面V形切口Fig.1.Bi-materialanti-planeV-notch（1）其中，是应力奇异性指数，是反平面V形切口在不同材料域内的位移特征角函数，是相应的位移组合系数，是截取的级数项数。双材料反平面V形切口位移-应变关系，可用位移角函数表示如下：（2）双材料反平面V形切口端部应力分量，可用位移角函数表示如下：（3）反平面应力平衡方程为：（4）将式（3）应力分量代入式（4）平衡方程，得（5）在和上的边界条件为（6a）（6b）在上的交界条件为：（6c）（6d）为了便于编程计算，引入无量纲变量，当时，，当时，则有，，，（7）将式（7）代入式（5），得（8）将式（7）代入式（6），得（9a）（9b）（9c）（9d）至此，双材料反平面V形切口端部区域应力奇异性指数的计算，转化为求解满足边界和交界条件式(9)时常微分方程组式(8)的特征值问题。3、微分求积法的基本理论图2双材料反平面V形切口的微分求积法计算模型Fig.2.Differentialquadraturemethodcalculationmodelofbi-materialantiplaneV-notch图2为双材料反平面V形切口的微分求积法计算模型，将两个楔形尖劈区间皆划分为段，考虑位移特征角函数在区间上可微，，个节点上函数的导数值可以用个节点上函数值的加权线性和近似表示，本文将函数用节点函数值进行拉格朗日(Lagrange)插值来描述，即(10)式中为拉格朗日多项式，其形式为(11)由式(10)对函数求一阶导数，得到(12)将式(12)在段区间上离散，从而进一步得到(13)将式(13)写成向量形式，得，(14)其中，显示表达式为(15)依次类推可得(16)这里，称为权系数矩阵。其中，权系数矩阵的各阶导数之间的关系为(17)将式（14）和（17）代入式（8），控制方程写成向量形式，得(18)其中，，，，，，，同理，边界和交界条件写成向量形式：(19a)(19b)(19c)(19d)这里，和。将边界和交界条件式（19）代入控制方程式（18），替换相应的0行和n行，得：(20)由于式（20）中含有项，导致非线性特征值问题，引入2个新向量和，得(21)其中，，联立式（20）和式（21）可得，(22)再令，，，则有，式(22)就简化为以下的矩阵特征方程组形式(23)至此，式(23)为个未知向量关于应力奇异性指数的一般广义代数特征方程组。对于广义代数特征值问题，本文将其转换成标准型特征值表达式。若存在逆矩阵，将式（23）变换为(24a)若存在逆矩阵，将式（23）变换为(24b)若和逆矩阵均不存在，但对于给定常数，存在，将式（23）变换为(24c)式（24）的代数特征值问题，可用数值方法求解，这里采用法求解该式，可得到双材料反平面V形切口端部应力奇异性指数和相应的特征向量。4、数值算例及问题讨论算例1.楔形尖劈张角β对反平面V形切口奇异性的影响如图1所示，作为双材料反平面V形切口的简单特例，当材料的剪切弹性模量之比时，双材料反平面V形切口就退化为各向同性单材料反平面V形切口，这里，从变化到，运用微分求积法，取区间的分段数4,8,12,16,20，且采用切比雪夫多项式的根作为区间上的布点方式时，本文计算值与文献[4]精确解同列于表1。表1计算结果表明：区间分段数越多本文方法计算精度越高，当计算值收敛于精确解，证明了微分求积法计算双材料反平面V形切口端部应力奇异性指数的有效性和计算精度。当两种材料的剪切模量比值时，，楔形尖劈的内张角从变化到，取区间分段数n=20，运用微分求积法计算，则双材料反平面V形切口奇性指数的计算值与文[4]和文[12]的计算结果同列于表2，从表2中可以看出，在所取的有效数字范围内，本文方法计算值能和文献[4]的精确解以及文献[12]有限元特征法获得的结果吻合。本文方法除了能计算出第1阶应力奇性指数以外，可以同时获得若干非奇异性高阶指数，表2共列出了本文方法计算的前3阶奇性指数，对于双材料反平面V形切口完整应力场来说，尽管和均大于1.0，不再呈现奇异性，但这些非奇异参数是描述切口尖端完整应力场不可或缺的量[15]。表1.当μ1/μ2=1时，奇性指数随尖劈角的变化Table1Variationofsingularityorderswithanglewhenμ1/μ2=1β精确解[4]微分求积法计算值（DQM）N=4N=8N=12N=16N=2010°-0.052631578-0.065660654-0.052635560-0.052631605-0.052631605-0.05263160530°-0.142857142-0.154645354-0.142860746-0.142857167-0.142857167-0-0.250000000-0.260314684-0.250003152-0.250000021-0.250000021-0.25000002190°-0.333333333-0.342501942-0.333336135-0.333333352-0.333333352-0.333333352120°-0.400000000-0.408251748-0.400002522-0.400000017-0.400000017-0.400000017150°-0.454545454-0.462047043-0.454547747-0.454545470-0.454545470-0.454545470180°-0.500000000-0.506876456-0.500002102-0.500000014-0.500000014-0.500000014表2.当μ1/μ2=10时，奇性指数随尖劈角的变化Table2Variationofsingularityorderswithanglewhenμ1/μ2=10Ref.[4]Ref.[12]PresentPresentPresent10°-0.00558-0.00558-0.0055801590.9884917801.9817769630°-0.01796-0.01796-0.0179615670.9485938541.7759955160°-0.04877-0.04877-0.0487688610.4999999581.0487687890°-0.13677-0.13677-0.1367776770.1367776210.99999994120°-0.28755-0.28755-0.2875504410.0445802320.95541968150°-0.41097-0.41096-0.4109684050.0165658400.75958238180°-0.50000-0.50000-0.5000000140.0000000000.49999995算例2.剪切弹性模量比值μ1/μ2对反平面V形切口奇异性的影响如图1所示，若，由文献[4]知，双材料反平面V形切口端部第1阶应力奇异性指数与两种材料楔形尖劈的内张角之和有关，而与其剪切弹性模量的比值无关，且其精确解为。本文取，从2.0变化到10.0，由微分求积法计算，双材料反平面V形切口应力奇异性指数与文献[4]精确解同列于表3，计算结果验证了文献[4]的结论。若，，时，双材料反平面V形切口的第1阶应力奇异性指数随着两种材料的剪切弹性模量的比值的变化规律如表4所示，由微分求积法计算，本文计算值与文献[4]精确解和文献[10]计算结果吻合，且随着比值增大，双材料反平面V形切口端部第1阶应力奇异性指数减小，即奇异性减弱。表3α=β=135°时，反平面V形切口应力奇异性指数随剪切弹性模量的比值μ1/μ2的变化Table3Variationofsingularityorderswithshearmodulusratioofbiomaterialsμ1/μ2whenα=β=135μ1/μ2精确解[4]微分求积法计算值（DQM）N=4N=8N=12N=162.0-0.333333333-0.333808531-0.333333359-0.333333352-0.3333333524.0-0.333333333-0.333808530-0.333333359-0.333333352-0.3333333526.0-0.333333333-0.333808530-0.333333359-0.333333352-0.3333333528.0-0.333333333-0.333808530-0.333333359-0.333333352-0.33333335210.0-0.333333333-0.333808530-0.333333359-0.333333352-0.333333352表4α=180°，β=90°时，反平面V形切口应力奇异性指数随剪切弹性模量的比值μ1/μ2的变化Table4Variationofsingularityorderswithshearmodulusratioofbiomaterialsμ1/μ2whenα=180°，β=90°μ1/μ2文献[4]文献[10]微分求积法计算值（DQM）切比雪夫布点方式N=4N=8N=12N=162.0-0.26772-0.267720468-0.270018139-0.267720657-0.267720493-0.2677204934.0-0.20483-0.204832758-0.208083351-0.204833115-0.204832787-0.2048327876.0-0.17223-0.172237321-0.176059070-0.172237812-0.172237354-0.1722373548.0-0.15147-0.151478016-0.155691333-0.151478610-0.151478048-00-0.13677-0.136777643-0.141281753-0.136778323-0.136777677-0.136777677式(1)中位移特征角函数是描述双材料反平面V形切口端部位移场和奇异应力场不可缺少的物理量，本文微分求积法既可计算出V形切口端部的应力奇异性指数，同时相应的应力特征角函数也一并解出，当，时，比值分别为4.0、10.0时，图3和图4是表4中双材料反平面V形切口第1阶主奇异特征值、所对应的位移特征角函数曲线图。由图3和图4可见，在异质材料粘结处位移场特征函数是连续的，而位移场一阶导函数在异质材料粘结处发生了突变，这一现象导致在双材料交界处应力分量相等而不相等，这个结果可能是导致界面破坏的原因之一。图3.双材料反平面V形切口位移特征角函数及其一阶导函数，μ1/μ2=4.0Fig.3.Displacementcharacteristicanglefunctionanditsfirstderivativefunctionforabi-materialV-notcheswhenμ1/μ2=4.0图4.双材料反平面V形切口位移特征角函数及其一阶导函数，μ1/μ2=10.0Fig.4.Displacementcharacteristicanglefunctionanditsfirstderivativefunctionforabi-materialV-notcheswhenμ1/μ2=10.05、结论基于双材料反平面V形切口端部附近区域位移场的Williams渐近展开式，给出了微分求积法分析切口端部应力奇异性计算的一种新途径，主要结论如下：基于线弹性力学理论，由本文微分求积法得到的双材料反平面V形切口端部应力奇性指数与精确解、已有文献的计算结果吻合，证明了微分求积法计算双材料反平面V形切口端部应力奇性指数的有效性。同时，本文方法具有计算量小、精度高的优势。微分求积法分析双材料反平面V形切口端部奇异应力场，可一次性解出双材料反平面V形切口端部所有的应力奇性指数及其相应的特征函数，计算结果发现：在不同材料粘结界面上位移场一阶导函数发生了突变，这就是多材料接头界面发生失效破坏的原因之一。在计算双材料反平面V形切口端部较远区域完整应力场时，V形切口端部非奇异应力奇性指数是不可或缺的物理量，本文方法同时获得了若干阶高阶非奇异应力奇性指数。基于本文微分求积法研究各向同性双材料反平面V形切口端部奇异应力场的基础，可将本文方法拓展到正交各向异性复合材料和压电材料反平面V形切口端部奇异应力场的分析研究中。参考文献[References]BogyDB.Twoedge-bondedelasticwedgesofdifferentmaterialsandwedgeanglesundersurfacetractions[J].TransactionofASME,JournalofAppliedMechanics,1971,38,377-386.DempseyJP,SinclairGB.Onthestresssingularitiesintheplaneelasticityofthecompositewedge[J].JournalofElasticity,1979,9,373-391.DempseyJP,SinclairGB.Onthestressbehaviouratthevertexofabi-materialwedge[J].JournalofElasticity,1981,11,317-327.MaCC,HourBL.Analysisofdissimilaranisotropicwedgessubjectedtoantiplanesheardeformation[J].InternationalJournalofSolidsandStructures,1989,25,1295-1309.TingTCT.Explicitsolutionandinvarianceofthesingularitiesataninterfacecrackinanisotropiccomposites[J].InternationalJournalofSolidsandStructures,1986,22,965-983.KargarnovinMH,ShahaniAR,FariborzSJ.Analysisofanisotropicfinitewedgeunderantiplanedeformation[J].InternationalJournalofSolidsandStructures,1997,34,113-128.Shahani,AR.Analysisofananisotropicfinitewedgeunderanti-planedeformation[J].JournalofElasticity,1999,56,17-32.平学成,陈梦成,谢基龙.复合材料尖劈和接头端部奇性场的反平面问题研究[J].固体力学学报,2005,26(2):193-198.(PingXC,ChenMC,XieJL.Singularstressstatesintipsofanisotropicmulti-materialwedgesandjunctionsubjectedtoantiplaneshear[J].ActaMechanicaSolidaSinica,2005,26(2):193-198.(inChinese))葛仁余,牛忠荣,程长征,胡宗军,杨智勇.反平面V形切口塑性应力奇异性分析[J].固体力学学报,2015,36(01):76-84.(GeRenyu,NiuZhongrong,ChengChangzheng.Analysisofplasticstresssingularitiesofanti-planeV-notchesinhardeningmaterials[J].ActaMechanicaSolidaSinica,2015,36(1):76-84.(inChinese))葛大丽.双材料反平面V形切口应力奇性指数的计算[J].合肥工业大学学报(自然科学版)