ANSYS粘弹体分析

1、 1. 1AXSYS中表征粘弹性属性问题 1 1.2Prony级数形式1 1. 3Maxwell 形式3 1.3建模与载荷条件5 1.3.1模型设计5 1.3.2 有限元建模5 1. 3. 3理论解析解计算式 6 1.4有限元数值解与结果比较 6 1. 4. 1 Planel83 Prony 级数方式 6 1.4.5算例结论 10 ANSYS中粘弹材质属性参数输入和分析 1. 1 ANSYS屮表征粘弹性属性问题 粘弹性材料的应力响应包括弹性部分和粘性部分，在载荷作用下弹性部分是即时响应的，而粘性部分需 要经过一段时间才能表现出来。一般的，应力函数是由积分形式给出的，在小应变理论下，各向同性的粘

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

3、axwell形式，一种是结构单元（如Planel83, Planel82等）所釆用的Prony级数形式。实 际上，这两种表示方式是一致的，只是具体数学表达式有一点点不同。 1.2 Prony级数形式 用Prony级数表示粘弹性属性的基本形式为： nGft） G （t ）二G+Gi exp.一 （0. 2） i 土s T 其中，c：和Gi是剪切模量，K和K】是体积模量，-i和-i是各Prony级数分量的松弛时间。再定义下面 相对模量 :=Gi ： Go (0. 4) (0. 5) 其中，Go, K。分别为粘弹性材质(固体推进剂)的瞬态模量，并定义式如下： Go G t = 0 =G . _-二

4、Gj(o. 6) i =1 Ko=K t=0 二K.亠二 Kx(o.7) i =1 GK 在ANSYS中，Prony级数的阶数口和0可以不必相同，当然其屮的松弛时间“和i也不必相同。 对于粘弹性问题，粘弹体的泊松比一般是取为时间的函数 - t。不过有时情况允许也可近似设为常 数，这时根据弹性常数关系就有： (0. 8) 3 1 -2J 其中，E t为松弛模量，由实验来确定。E t ,G t ,K t的相应系数比相同。 这样就可以将G t和K t统一于E t形式。若我们将松弛模量表示为Prony级数形式，即: n f t ) E t、Ei exp 1 1 (0. 9) 1匸丿 i =1 G K

5、于是，G t和K t中有， n .J = i = i y 寸fko类似于G。、 Ko ,我们也同样定 义瞬态松弛模量E。： EE t =0i=E：: Ex (0. 10) i甘 这样，由(0.8)可得 Go亠 2 1J (0. 11) Eo Ko = 3 1 -2 J 要注意的是，ANSYS中对Prony级数的支持项数不能超过6项，即n乞6。这确实是一个遗憾。 另外， The viscoelasticity in put for SHELL181 , PLANE182 , PLANE183 , SOLID185 , SOLID186 , SOLID187 , S0LSH190 , SHELL2

6、08 , and SHELL209 consists of elasticity properties and relaxation properties The underlying elasticity is specified by either the MP comma nd (for hypoelasticity) or by the TB , HYPER comma nd (for hyperelasticity) Use the TB , PRONY or TB , SHIFT comma nds to i nput the relaxation properties 可见，此时

7、除了由Prony级数形式附加粘弹性，还需输入弹性”属性。这里我对hypoelasticity T ，具 体也说不上来。在ANSYS帮助文档里有这样一段： !Small Strain Viscoelasticity mp,ex,1,20. 0E5 !elastic properties mp, nu xy, 1, 0.3 tb, pr ony , 1, , 2, shear tbdata, 1, 0. 5, 2. 0, 0. 25, 4. 0 !defi ne viscosity parameters (shear) tb, pro ny, 1, 2, bulk !defi ne viscosi

8、ty parameters (bulk) tbdata, 1, 0. 5, 2. 0, 0. 25, 4. 0 !Large Strain Viscoelasticity tb, hyper, 1, , m oon tbdata, 1, 38. 462E4, , 1. 2E6 !elastic properties tb, pr on y, 1, , 1, shear tbdata, 1, 0. 5, 2. 0 !defi ne viscosity parameters tb, pro ny, 1, 1, bulk tbdata, 1, 0. 5, 2. 0 !defi ne viscosit

9、y parameters 1.3 Maxwell 形式 For the viscoelastic elements VISC088and VISC089the material properties are expressed in integral form using the kernel function of the generalized Maxwell eleme nts as: GgG Gexp nG (0. 12) nK K =K .二 Ki exp i _1 Cz =G Go-G (0. 13) Di 二Ki/ Ko -K：; 其中F为折算时间，由于不考虑温度载荷，方程中的折

10、算时间就是实际时间，即情形的。 二仁类同Prony级数 E = reduced or peudo time G( E ) = shear relaxation kernel function K( E ) = bulk relaxation kernel function nG = nu mber of Maxwell eleme ntsusedtoapproximate the shear relaxati on kern el (input con sta nt 50) nK = numb erof Maxwell eleme nts used to approximate thebulk

11、relaxatio n kern el (input con sta nt 71) Ci =con sta nts associated withthein sta ntan eous resp onse for shear behavior(in putcon sta nts 51 -60)Di = constants associated with the instantaneous responsefor bulk behavior (input constants 76 -85)Go =initial shear modulus(in put con sta nt 46) G ： ：

12、= final shear modulus (in put con sta nt 47) Ko = in itial bulk modulus (in put con sta nt 48) K ： ： = final bulk modulus (input constant 49) G i = con sta nts associated with a discrete relaxati on spectrum in shear (in put con sta nt 61-70) K i = con sta nts associated with a discrete relaxati on

13、spectrum in bulk (in put con sta nt 86-95) 同Prony技术情形一样的： 由试验数据拟合得到(0.12)； 由(o.i2)即可确定：级数项数a, nK ；k和g的初始值和稳态值：Kc, K：:和G,g：；时间松弛系数*、； 再分别根据(0.13)计算得到参数G, Djo 将上面计算所得值分别填入Maxwell材质属性表即可。 Here, Go and K o are, respectively, the shear and bulk moduli at the fast load limit (ie. the in sta ntan eous modu

14、li), and G. 一 and i are the moduli at the slow limit The elasticity parameters in put corresp ond to those of the fast load limit. In itialize the con sta nt table with TB , EVISC You can define up to 95 con sta nts (C1-C95) with TBDATA comma nds (6 per comma nd)： 1.3建模与载荷条件1.3.1模型设计 如图3.1-1所示，一个圆孔形

15、的药柱，内径为 6外径为b,弹性钢壳体厚度为ho药柱内表面受均布压 强载荷作用。另外，我们假设： a.药柱外表面与壳体是直接粘接在一起的，忽略绝热层等材料的厚度；b、 该圆孔型药柱足够长，可以简化为平面应变问题来处理；c、推进剂泊松比为常数。其中，相关物性参数如 下： 1、壳体：E=196.5 Gpa,V二0.29,并认为v为常数处理； 2、推进剂:V =0. 495,松弛模量E(t)用Prony级数表示为： _t_t_t E(t)二0.7058860. 168169 e1H；0. 098714 e3013j 07 1. 930384 e30b 307 (MPa) (0. 14) 3、几何参数

16、：a= 100mm, b=177mm, h=3mm 。 4、阶跃压力载荷： P t 二Fo Id(0. 其中,P。为稳态压强值,取/F=6.3238MPa, n =20o我们还可以将所得结果于其加以对比。 1.3. 2有限元建模 图3. 1-1所示的模型的几何形状以及所受载荷条件和边界条件具有明显的对称性，出于方便建模和适当 计算量的考虑，我们取其1/4来建模分析。根据前面药柱外表面与壳体内表面是直接粘接在一起的假设，可 以将这两个面位移耦合。其有限元模型和网格划分如图3. 1-2所示。图中还标岀模型位移边界约束，内表面 压强载荷以及药柱与壳体粘结面上节点耦合约束。 这里我们统一选用Plane

17、l83单元来划分网格的，共200个单元，718个节点：其 中，壳体部分有单元50个，节点205个；药柱部分有单元150个，节点513个。 AN 图 3.1-2 求解时间定为iso 下面就不同加载方式，分别计算分析。我们将就节点A (如图3. 1-2中所示)取其径向应力/应变和周向 应力/应变加以比较。 1.3. 3理论解析解计算式 限于篇幅, A2 这里不准备列写出理论解析解的详细公式推导过程，根据文献，这里直接给岀最终计算表达式。 径向应力: 径向应变: 1 2vv 1 1 -2v、2 L -2v ) (2)周向应力: Jr, -A2 2 a 1+ 周向应变: 2v 1 v 1 /二上 Po

18、 i +(i -2v Tfrj *2 h Ec b “2 mk-c2) 其中，a. b为药柱的内径和外径，h为壳体厚度，怎二b/a, E,为壳体的弹性模量, 1 F s E s。即Ft可由Et作一定变换 F t是推进剂蠕变柔量，与松弛模量 E t有如下关系: 得到。 1.4有限元数值解与结果比较 1. 4. 1 Planel83 Prony级数方式参数变换 另外，根据前面Prony级数表示方式，经换算得到相关各系数为: Eo =2. 903153MPa, v =0. 495 常数;根据(0.8)式， :二；=30130. 7,二0.0579 I 二 2 =3013. 07, : ? 一鳥二0.

19、 0340 G 二.K =301. 307,：3“ 二：-K =0. 6649 er i d. No dl Huftiber 1 L LL MQIIAlE Aval 1 abl# 参数输入情况分别如下图所示: Lmeiar Isstr pic Prijjiy : ShfVA* BAEP0D*4 Prcny s Vi-lrnfttric trt H-a. tei* 5 V N tidel Bunbeif 2. Linear Isot ropxc Propert ies for Eater X iLin 电世 I 弓otropit MaleriAl Pr&ptrlies *r Material Hsb r t 9 T1 EX FEXT 2 9032E006 Kdd Tenperature jDelete Temperalwre 1 orCancel &r*ph Help 1 NtoriModels Material Model Ntjtiber I LirwH Isotropic Frony: Show Rypon Prony： Voliwetric Res 5 Nateril Model Num