ANSYS粘弹体分析

1、ANSYS中粘弹材质属性参数输入和分析1.1 ANSYS中表征粘弹性属性问题11.2 Prony级数形式11.3 Maxwell形式31.4 建模与载荷条件51.4.1 模型设计51.4.2 有限元建模51.4.3 理论解析解计算式61.5 有限元数值解与结果比较61.5.1 Plane183,Prony级数方式61.5.5 算例结论10ANSYS中粘弹材质属性参数输入和分析1.1ANSYS中表征粘弹性属性问题粘弹性材料的应力响应包括弹性部分和粘性部分，在载荷作用下弹性部分是即时响应的，而粘性部分需要经过一段时间才能表现出来。一般的，应力函数是由积分形式给出的，在小应变理论下，各向同性的粘弹性

2、本构方程可以写成如下形式：(0.1)042G(t-喷其中c=Cauchy应力G（t）=为剪切松弛核函数K（t）=为体积松弛核函数6=为应变偏量部分（剪切变形）=为应变体积部分（体积变形）=当前时间七=过去时间=为单位张量。该式是根据松弛条件本构方程（0.1）,通过将一点的应变分解为应变球张量（体积变形）和应变斜张量（剪切变形）两部分，推导而得的。这里不再敖述，可参考相关文献等。ANSYS中描述粘弹性积分核函数G（t.口K（t谬数表示方式主要有两种，一种是广义Maxwell单元（VISCO88和VISCO89）所采用的Maxwell形式，一种是结构单元（如Plane183,Plane182等）所

3、采用的Prony级数形式。实际上，这两种表示方式是一致的，只是具体数学表达式有一点点不同。1.2 Prony级数形式用Prony级数表示粘弹性属性的基本形式为：(0.2)nGGt=G.Giexpi1nKtKt=K.KiexpI穴i2,i(0.3)其中，G毛和Gi是剪切模量，Ka和Ki是体积模量，7G和是各Prony级数分量的松弛时间。再定义下面相对模量(0.4)'2i-KiK0(0.5)其中，Go,Ko分别为粘弹性材质(固体推进剂)的瞬态模量，并定义式如下:nGG0=Gt=0=GJ-,Gii1(0.6)nKK0=Kt=0)=K二八Kii1(0.7)在ANSYS中，Prony级数的阶数n

4、G和nK可以不必相同，当然其中的松弛时间和。和却K也不必相同。对于粘弹性问题，粘弹体的泊松比一般是取为时间的函数N=N(t卜不过有时情况允许也可近似设为常数，这时根据弹性常数关系就有:(0.8)EtKt)31-2其中，E(t)为松弛模量，由实验来确定。E(t»G(t),K(t)的相应系数比相同。这样就可以将G(t)和K(t)统一于E(t)形式。若我们将松弛模量表示为Prony级数形式，即:E(t产Eg+ZEiexp-(0.9)是，G(t)和K(t)中有，n=nG=nK,a=TiK,必=5g=5k。类似于G。、K0,我们也同样定义瞬态松弛模量E0:(0.10)nGE。=Et=0)=E：

5、'.旦i1这样，由(0.8)可得Eo一21(0.11)Eo一312J要注意白是，ANSYS中对Prony级数的支持项数不能超过6项，即nW6。这确实是一个遗憾。另外，TheviscoelasticityinputforSHELL181,PLANE182,PLANE183,SOLID185,SOLID186,SOLID187,SOLSH190,SHELL208,andSHELL209consistsofelasticitypropertiesandrelaxationproperties.TheunderlyingelasticityisspecifiedbyeithertheMPcom

6、mand(forhypoelasticity)orbytheTB,HYPERcommand(forhyperelasticity).UsetheTB,PRONYorTB,SHIFTcommandstoinputtherelaxationproperties.可见，此时除了由Prony级数形式附加粘弹性，还需输入“弹性”属性。这里我对hypoelasticity不了解，具体也说不上来。在ANSYS帮助文档里有这样一段：!SmallStrainViscoelasticitymp,ex,1,20.0E5!elasticpropertiesmp,nuxy,1,0.3tb,prony,1,2,shear

7、!defineviscosityparameters(shear)tbdata,1,0.5,2.0,0.25,4.0tb,prony,1,2,bulk!defineviscosityparameters(bulk)tbdata,1,0.5,2.0,0.25,4.0!LargeStrainViscoelasticitytb,hyper,1,moon!elasticpropertiestbdata,1,38.462E4,1.2E-6tb,prony,1,1,shear!defineviscosityparameterstbdata,1,0.5,2.0tb,prony,1,1,bulk!define

8、viscosityparameterstbdata,1,0.5,2.01.3 Maxwell形式FortheviscoelasticelementsVISCO88andVISCO89thematerialpropertiesareexpressedinintegralformusingthekernelfunctionofthegeneralizedMaxwellelementsas:nGG(亡户G"Gexp.一，工"J(0.12)nKf七)K(5=%+£Kexp会:i工,/tjfjZG=GiG°-G二(0.13)，=t,类同Prony级数D=KiKo-

9、K:;其中己为折算时间，由于不考虑温度载荷，方程中的折算时间就是实际时间，即情形的。E=reducedorpeudotimeG(E)=shearrelaxationkernelfunctionK(E)=bulkrelaxationkernelfunctionnG=numberofMaxwellelementsusedtoapproximatetheshearrelaxationkernel(inputconstant50)nK=numberofMaxwellelementsusedtoapproximatethebulkrelaxationkernel(inputconstant71)Ci=c

10、onstantsassociatedwiththeinstantaneousresponseforshearbehavior(inputconstants51-60)Di=constantsassociatedwiththeinstantaneousresponseforbulkbehavior(inputconstants76-85)G0=initialshearmodulus(inputconstant46)G二二finalshearmodulus(inputconstant47)K0=initialbulkmodulus(inputconstant48)K：=finalbulkmodul

11、us(inputconstant49)Gi=constantsassociatedwithadiscreterelaxationspectruminshear(inputconstant61-70)K'i=constantsassociatedwithadiscreterelaxationspectruminbulk(inputconstant86-95)同Prony技术情形一样的：由试验数据拟合得到(0.12);由(0.12)即可确定：级数项数nG,nK；K和G的初始值和稳态值：K0,K0c和G0,G*;时间松弛系数、九K；再分别根据(0.13)计算得到参数Ci,DiO将上面计算所得

12、值分别填入Maxwell材质属性表即可。Here,G0andK0are,respectively,theshearandbulkmoduliatthefastloadlimit(i.e.theinstantaneousmoduli),andG二一andK.arethemoduliattheslowlimit.Theelasticityparametersinputcorrespondtothoseofthefastloadlimit.InitializetheconstanttablewithTB,EVISC.Youcandefineupto95constants(C1-C95)withTBD

13、ATAcommands(6percommand):1.3 建模与载荷条件1.3.1 模型设计如图3.1-1所示，一个圆孔形的药柱，内径为a,外径为b,弹性钢壳体厚度为ho药柱内表面受均布压强载荷作用。另外，我们假设：a、药柱外表面与壳体是直接粘接在一起的，忽略绝热层等材料的厚度；b、该圆孔型药柱足够长，可以简化为平面应变问题来处理；c、推进剂泊松比为常数。其中，相关物性参数如下：1、壳体：E=196.5Gpa,v=0.29,并认为v为常数处理；2、推进剂：v=0.495,松弛模量E(t)用Prony级数表示为：_t_t_tE(t)=0.7058860.168169e70.098714e3013

14、.071.930384e301.307(MPa)(0.14)3、几何参数：a=100mm,b=177mm,h=3mm。图3.1-14、阶跃压力载荷:Pt=811(0.15)其中，Po为稳态压强值，取R=6.3238MPa,n=20。我们还可以将所得结果于其加以对比。1.3.2 有限元建模图3.1-1所示的模型的几何形状以及所受载荷条件和边界条件具有明显的对称性，出于方便建模和适当计算量的考虑，我们取其1/4来建模分析。根据前面药柱外表面与壳体内表面是直接粘接在一起的假设，可以将这两个面位移耦合。其有限元模型和网格划分如图3.1-2所示。图中还标出模型位移边界约束，内表面压强载荷以及药柱与壳体粘

15、结面上节点耦合约束。这里我们统一选用Plane183单元来划分网格的，共200个单元，718个节点：其中，壳体部分有单元50个，节点205个；药柱部分有单元150个，节点513个。图3.1-2求解时间定为1s。下面就不同加载方式，分别计算分析。我们将就节点A(如图3.1-2中所示)取其径向应力/应变和周向应力/应变加以比较。1.3.3 理论解析解计算式限于篇幅，这里不准备列写出理论解析解的详细公式推导过程，根据文献，这里直接给出最终计算表达径向应力:;二r,t二U1W'2-1_.r2b2.12.r2I.1-2*E±hEc径向应变:九2(1+P0(2)周向应力:三r.tPo2-

16、A-2-1S1-2r2-2-1r2bELhEc周向应变:"-ItfP0£(1+v此a2L/.|(1-2v)+P2-1_b1hEc其中,a、b为药柱的内径和外径，h为壳体厚度，？=b/a,Ec为壳体的弹性模量,(1c2)1一.F(t心推进剂蠕变柔量，与松弛模量E(t)有如下关系：F(sF(s尸/。即F(t)可由E(t)作一定变换得到。1.4 有限元数值解与结果比较1.4.1 Plane183、Prony级数方式参数变换另外，根据前面Prony级数表示方式，经换算得到相关各系数为:巳=2.903153MPa,v=0.495常数;1G=；=30130.7,:1G=二；=0.057

17、9.22=K=3013.07,寸-:K=0.0340,G=.；=301.307,：-3g=:；=0.6649根据（0.8）式,参数输入情况分别如下图所示:IodelBehaviororiMb.t«nblNadalsDa-tini-dlu.MqiiaIeAval1Abla-literkalNodwl5ber1FilLinearIsotropxcPropertiesforEater.X|LtnearIs-strpicPrijjiy:Shfl'A"&曲印5£一Preny：VolumetricRespoiLin电旺IstklrftpieMalerid.Pr

18、&ptrlies£*rMaterial耳画后总工I-tertHaEWiidLNtidelBumbeif2.TlEXFEXT29032E>006Add.TenperaturejDeleteTefliperalwreOKCancelHelpLIattrialMedelVufib«rPronyTablelx卷LinearIsetrcrpich*4¥：SLt«rR”Qn$ia号Prony：VolurictfLC现上logNaterialModelffunber2TlTemperalwr电al11PC131a2D.034002隹soial006649

19、3t3SOL31AddBowDeleteRxwAdidTewRer刊tor电口rl电teTmuer3tMeJPMHYTitleforH&teraaltfunber1-PronyVicoelasticShearResp<?ELse计算结果比较计算值与理论值比较，结果如下面几个图中所示以圆通内表面点A的应力应变为参考来考察，将ANSYS的：周向应变%才0)0061*，V*耙*/理也解掰解+ansH#值解片向应力斗网尸蓟-理徽解析初+ANS彩效值解周向应力a+AMSM值解7口.日0.911.4.2 Maxwell形式参数变换与输入根据前面Prony级数表示方式，由(0.12)即可确定：

20、级数项数nG=nK=3；Ko=9.6772e007,K笛=2.3530e007和Go=9.7095e+005,Gk2.3608e+005;时间松弛系数iG=%K=30130.73013.07301.307；再分别根据(0.13)计算得到参数Ci=Di=0.07650.04490.8785o因为我们这里是假设泊松比是常数的；否则，拟合所得K(t网G(t)形式就没有(0.8)那样简洁的关系，当然E(t)G(t),K(t)的相应系数比也就不相同了。而且nG,nK亦可以取为不同，这都不影响上面参数的换算的。参数输入情况分别如下图所示。其中无需对EX等参数输入。ViscDclasticityforlat

21、erialBimbcr1Visco-elas-ticityforMaterialNumber1Cl'CS12345QjcTFC6-CWp_LpVC11-C15FVVC16-C20F|oxxC21-C25rb卜C26-C30FrkC31-C3SFkVVxC3&-C40FTTkC41-C45FkVVVC4&-C50忸TM5E*005k3M6E-H0O5|9BlTSEtOOTJ2.353E-KKTFC51-C55|0.(776536b0M926|0.87B54XpC5B-C80Fb丁C61-C65|30131bo131|301.31bC66-C70FTx"FC71=C75rFkV期CT6-C80|o076536044926|o37854VVC61-C65丁|oTkkC86-C90boi3i13cli31|30131xxC91-C95rPr>T计算结果比较同样以圆通内表面点A的应力应变为参考来考察，将ANSYS计算值与理论值比较，结果如下面几个图中所示的：0.0811十尸千十十中中+节*Illi1i11il-1/f11理希第祈解/1+时喈年其向解r：i:冏向应变电正外0.040.0200.1D.2D.3Q.40.50.6D70月0.910rT1:i1：1