ANSYS粘弹体分析

1、ANSYS中粘弹材质属性参数输入和分析 11.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.1 ANSYS中表征粘弹性属性问题粘弹性材料的应力响应包括弹性部分和粘性部分， 在载荷作用下弹性部分是即时响应的， 而粘性部分需 要经过一段时间才能表现出来。 一般的，应力函数是由积分形式给出的， 在小应变

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

3、是结构单元（如 Plane183, Plane182等）所采用的Prony级数形式。实际上，这两种表示方式是一致的，只是具体数学表达式有一点点不同。1.2 Prony 级数形式用Prony级数表示粘弹性属性的基本形式为：nGr t ）G （t 尸G.+2 Gi exp - ；（0.2）i 二LT(0.3)nKr t)K (t 尸Kg+£ Ki exp . f -13其中，G毛和Gi是剪切模量,K3c和Ki是体积模量，7G和qK是各Prony级数分量的松弛时间。再定义下面相对模量(0.4)-：SG =Gi；Go-：i = Ki . Ko(0.5)其中，G0, K0分别为粘弹性材质(固体

4、推进剂)的瞬态模量，并定义式如下:nGG0 =G t = 0 =G，-二 Gi(0.6)i =1 nKK0 = K t =0 = K-=： - Ki(0.7)1 4在ANSYS中，Prony级数的阶数nG和nK可以不必相同，当然其中的松弛时间丁和£，也不必相同。对于粘弹性问题，粘弹体的泊松比一般是取为时间的函数P = k(t)。不过有时情况允许也可近似设为常数，这时根据弹性常数关系就有：E tG t )2 1 口(0.8)E t K t )3 1 -2其中，E(t)为松弛模量，由实验来确定。E(t ),G(t ),K(t )的相应系数比相同。这样就可以将G (t )和K (t )统一

5、于E (t )形式。若我们将松弛模量表示为Prony级数形式，即：n'' t、E (t )=E3O+X Ei exp (0.9)日i J于是，G(t)和K(t)中有， n=nG=nK,%=%G=TiK, 5 =c(iG=砰。类似于G0、K0,我们也同样定义瞬态松弛模量E0: nG E0 = E t = 0 =E=% Ei(0.10)1 1这样，由(0.8)可得GoE02 1(0.11)KoE03 1 -2 J要注意白是，ANSYS中对Prony级数的支持项数不能超过 6项，即n <6 o这确实是一个遗憾。另外,The viscoelasticity input for S

6、HELL181 , PLANE182 , PLANE183 , SOLID185 , SOLID186 , SOLID187 , SOLSH190 , SHELL208 , and SHELL209 consists of elasticity properties and relaxation properties . The underlying elasticity is specified by either the MP command (for hypoelasticity) or by the TB ,HYPER command (for hyperelasticity). Us

7、e the TB ,PRONY or TB ,SHIFT commands to input the relaxation properties.可见，此时除了由 Prony级数形式附加粘弹性，还需输入“弹性”属性。这里我对hypoelasticity不了解，具体也说不上来。在 ANSYS帮助文档里有这样一段：!Small Strain Viscoelasticitymp,ex,1,20.0E5 !elastic propertiesmp,nuxy,1,0.3tb,prony,1,2,shear!define viscosity parameters (shear)tbdata,1,0.5,2

8、.0,0.25,4.0tb,prony,1,2,bulk!define viscosity parameters (bulk)tbdata,1,0.5,2.0,0.25,4.0!Large Strain Viscoelasticitytb,hyper,1,moon!elastic propertiestbdata,1,38.462E4,1.2E-6tb,prony,1,1,shear!define viscosity parameterstbdata,1,0.5,2.0tb,prony,1,1,bulk!define viscosity parameterstbdata,1,0.5,2.01.

9、3 Maxwell 形式For the viscoelastic elements VISCO88and VISCO89the material properties are expressed in integral form using the kernel function of the generalized Maxwell elements as:nGB( t )G ( _ )=G3c+£ Gi exp .(0.12)nKf七)K (V)=K0c+£ Ki exp -1- 1九 JCi =GiGo-G:;(0.13)D =K Ko -K:;Prony级数其中己为折

10、算时间，由于不考虑温度载荷，方程中的折算时间就是实际时间，即 U=t,类同 情形的。E = reduced or pseudo timeG( E ) = shear relaxation kernel functionK( E ) = bulk relaxation kernel functionnG = number of Maxwell elements used to approximate the shear relaxation kernel (input constant 50)nK = number of Maxwell elements used to approximate

11、the bulk relaxation kernel (input constant 71)Ci = constants associated with the instantaneous response for shear behavior (input constants 51 -60)Di = constants associated with the instantaneous response for bulk behavior (input constants 76 -85)Go = initial shear modulus (input constant 46)G 二二fin

12、al shear modulus (input constant 47)Ko = initial bulk modulus (input constant 48)K ： = final bulk modulus (input constant 49)Gi = constants associated with a discrete relaxation spectrum in shear (input constant 61-70)Ki = constants associated with a discrete relaxation spectrum in bulk (input const

13、ant 86-95)同Prony技术情形一样的：由试验数据拟合得到(0.12);由(0.12)即可确定：级数项数nG,nK ； K和G的初始值和稳态值：K0, K0c和G0,G*;时间松弛系数再分别根据(0.13)计算得到参数Ci, Dio将上面计算所得值分别填入Maxwell材质属性表即可。Here, Go and K 0 are, respectively, the shear and bulk moduli at the fast load limit (i.e. the instantaneous moduli), andG二一 and K.are the moduli at the

14、slow limit. The elasticity parameters input correspond to those of the fast load limit.Initialize the constant table with TB ,EVISC. You can define up to 95 constants (C1-C95) with TBDATA commands (6 per command):1.3建模与载荷条件1.3.1 模型设计如图3.1-1所示，一个圆孔形的药柱，内径为 a,外径为b,弹性钢壳体厚度为 ho药柱内表面受均布压 强载荷作用。另外，我们假设：a、

15、药柱外表面与壳体是直接粘接在一起的，忽略绝热层等材料的厚度；b、该圆孔型药柱足够长，可以简化为平面应变问题来处理；c、推进剂泊松比为常数。其中，相关物性参数如下：1、壳体： E=196.5 Gpa, v=0.29,并认为v为常数处理；2、推进剂：v=0.495,松弛模量E(t)用Prony级数表示为： _t_t_tE(t) =0.705886 +0.168169 e30130.7 +0.098714 e3013.07 +1.930384 e301.307 (MPa) (0.14)3、几何参数：a=100mm, b=177mm, h=3mm 。图 3.1-14、阶跃压力载荷:P t 1 -e,t

16、(0.15)其中，R为稳态压强值，取 R =6.3238MPa , n = 20。我们还可以将所得结果于其加以对比。1.3.2 有限元建模图3.1-1所示的模型的几何形状以及所受载荷条件和边界条件具有明显的对称性，出于方便建模和适当计算量的考虑，我们取其 1/4来建模分析。根据前面药柱外表面与壳体内表面是直接粘接在一起的假设，可 以将这两个面位移耦合。其有限元模型和网格划分如图3.1-2所示。图中还标出模型位移边界约束，内表面压强载荷以及药柱与壳体粘结面上节点耦合约束。这里我们统一选用Plane183单元来划分网格的，共 200个单元，718个节点：其中，壳体部分有单元 50个，节点205个；

17、药柱部分有单元 150个，节点513个。图 3.1-2求解时间定为1s。下面就不同加载方式，分别计算分析。我们将就节点A (如图3.1-2中所示)取其径向应力/应变和周向应力/应变加以比较。1.3.3 理论解析解计算式限于篇幅，这里不准备列写出理论解析解的详细公式推导过程，根据文献，这里直接给出最终计算表达(1)径向应力:径向应变:T 7 )Po21 11 -2v 1-b7 1 F t ' 21V 1-2v -当r J ， 工 一1_r(2)周向应力:Po一 r2r h Ec(4)周向应变:7 = _J-2v 1 v 弓Po11 -2v2 r2-1F t . . J ； 1-2v 为A

18、/-1_r h Ec其中，a、b为药柱的内径和外径，h为壳体厚度，？ = b/a,Ec为壳体的弹性模量，2 一 , 2(1 - ):2(-" 1)1 ' ''2?二-,-二-(1 fc )1 ' (1 -2'.) ,1 ' (1 -2'.) , 1 f.1 一.F t心推进剂蠕变柔量，与松弛模量E t有如下关系：F(s(s)=-2。即F t )可由E t作一定变换s得到。1.4 有限元数值解与结果比较1.4.1 Plane183、Prony 级数方式 参数变换另外，根据前面 Prony级数表示方式，经换算得到相关各系数为:Eo

19、=2.903153MPa, v=0.495 常数;1G = ； =30130.7, : 1G =二：=0.0579GKGK,2 = 2 =3013.07,、'12 -：'2 =0.0340,G =.3 =301.307,1 3G =:；=0.6649根据（0.8）式,参数输入情况分别如下图所示:Iodel Behavioror iMb.t«n bl Nadals Da-tini-dl u. MqiiaIe Aval 1 Abla-literk al Nodwl 5ber 1Fil Linear Isot ropxc Propert ies for Eater. . .

20、 X |Lt near Is-strpicPrijjiy: Shfl'A" &曲印5£一Preny： Volumetric Respoi-tert HaEWiidL N tidel Bumbeif 2.at trial Medel Vufib«r卷 Linear Isetrcrpic h*4¥： SLt«r R”Qn$ia号 Prony： VolurictfLC 现上log Nat eri al Model ffunber 2Lin电旺 Istklrftpie Malerid. Pr&ptrlies £*r M

21、aterial 耳画后总工 ITlEXFEXT2 9032E>006Add. TenperaturejDelete TefliperalwreOKCancelHelp11LIProny TableIXTlTemper alwr 电 al11PC131a20.034002隹so ia l00 66493t330L 31Add Bow Delete RxwA did TewRer 刊 tor 电口 rl 电 te Tmuer 3tMeJPMHY Title for H&tera al tfunber 1 - Prony Vicoelastic Shear Resp<?ELse0

22、cel Hep计算结果比较以圆通内表面点 A的应力应变为参考来考察，将ANSYS计算值与理论值比较，结果如下面几个图中所示的:周向应变%才0) 0 061*，V* 耙*/理也解掰解+ ansH# 值解CD.5A.3O 丁a.aa片向应力斗网尸蓟-理徽解析初+ ANS彩效值解周向应力a+ AMSM值解7 口.日 0.911.4.2 Maxwell 形式参数变换与输入根据前面Prony级数表示方式，由(0.12)即可确定：级数项数nG =nK =3；K0 =9.6772e007, KaO = 2.3530e007 和 G0 =9.7095e+005, G 2.3608e+005 ;时间松弛系数 7

23、； ="K =30130.7 3013.07 301.307；再分别根据(0.13)计算得到参数Ci = Di =0.0765 0.0449 0.8785 o因为我们这里是假设泊松比是常数的；否则，拟合所得K(tDG(t )形式就没有(0.8)那样简洁的关系，当然E (t )G (t )K (t )的相应系数比也就不相同了。而且nG ,nK亦可以取为不同，这都不影响上面参数的换算的。参数输入情况分别如下图所示。其中无需对EX等参数输入。ViscDclasticity for laterial Bimbcr 1Visco-elas-ticity for Material Number

24、1Cl'CS12345QjcTFC6-CWp_LpVC11-C15FVVC16-C20F|oxxC21-C25rb卜C26-C30FrkC31-C3SFkVVxC3&-C40FTTkC41-C45FkVVVC4&-C50忸 TM5E*005k3M6E-H0O5|9 BlTSEtOOTJ2.353E-KKTFC51-C55|0. (776536b0M926|0.87B54XpC5B-C80Fb丁C61-C65|30131bo 13 1|301.31bC66-C70FTx"FC71=C75rFkV期CT6-C80|o 076536044926|o 37854VVC61-C65丁|oTkkC86-C90boi3i13 cli3 1|301 31xxC91-C95rPr>T计算结果比较同样以圆通内表面点 A的应力应变为参考来考察，将 ANSYS计算值与理论值比较，结果如下面几个图中 所示的：0.0811十尸千十十中中+节*Illi 1 i 1 1 il - 1/f11理希第祈解/1+时喈年其向解r ：i :冏向应变电正外0.040.