华科大压电材料.doc

一、Y切石英板的静态分析1.1 问题描述计算对象为一块石英板，其尺寸为：长度方向（x轴）：L=0.5m厚度方向（y轴）：H=2h=0.01m宽度方向（z轴）：W=0.05m板厚度方向两个平面有电极，大小为V/2=1000v。用ANSYS计算自由和夹支边界两种情况下的电极所在平面单位面积上的电荷e。该石英板的材料常数为：密度=2650kg/m3弹性常数矩阵：压电常数矩阵：介电常数矩阵：1.2 理论解边值问题为：电极所在平面单位面积上电荷的理论解为：（1）自由表面（2）夹支表面1.3 ANSYS求解(1) 单元选择根据该石英板的材料属性，可以选择solid226单元进行求解。Solid226单元为耦合场20节点六面体单元，KEYOPT(1)=1001即可激活压电自由度、位移和电压。输入“/prep7”进入前处理。APDL：et,1,solid226,1001ET命令的功能为从单元库中指定一个单元类型。其使用格式为：ET, ITYPE, Ename, KOP1, KOP2, KOP3, KOP4, KOP5, KOP6, INOPR(2) 介电系数矩阵由于该石英板介电常数为不是各项同性的，所以需要用TB，DPER和TBDATA命令指定各向异性介电系数矩阵中的各项，ANSYS中需要输入的是相对介电常数，因此需要用介电常数除以真空中介电常数0=。TBOPT=0则指定恒应变介电系数矩阵S； TBOPT=1则指定恒应力介电系数矩阵T，后者将被转换为恒应变介电系数矩阵S。值得指出的是，solid5单元只能用MP命令制定PERX，PERY和PERZ三项，所以本次计算不采用。ANSYS中输入相对介电系数矩阵的各项时也需要注意顺序，其顺序如下：从而该石英的介电系数矩阵可以按如下输入：tb,DPER,1tbdata,1,4.4285,4.4974,4.5642tbdata,5,0.0971TB和TBDATA指令用于输入矩阵中的各项。其中DPER用于指定该矩阵为各向异性相对介电系数矩阵，另外ANEL为各项异性弹性系数矩阵，PIEZ为压电系数矩阵。与ANSYS中介电系数矩阵的输入类似，弹性系数矩阵和压电系数矩阵的输入也要遵从ANSYS中的规则。(3) 弹性系数矩阵弹性系数矩阵主要包括刚度系数矩阵c和柔度系数矩阵s。用户可以用TBOPT在这两种矩阵中选择。ANSYS中刚度系数矩阵c中常数的顺序为：从而该石英的弹性系数矩阵可以按如下输入：tb,ANEL,1tbdata,1,86.74E9,-8.25E9,27.15E9,-3.66E9 ! c11,c21,c31,c61,c41tbdata,7,129.77E9,-7.42E9,5.7E9 ! c22,c32,c62,c42tbdata,12,102.83E9,9.92E9 ! c33,c63,c43 tbdata,16,29.01E9,2.53E9 ! c66,c46,c56tbdata,19,38.61E9 ! c44tbdata,21,68.81E9 ! c55(4) 压电系数矩阵压电矩阵可以定义压电应力矩阵e或压电应变矩阵d。压电应力矩阵e与各项异性弹性刚度矩阵c相联系，而压电应变矩阵d与弹性柔度矩阵s相联系。ANSYS中压电应力矩阵e中常数的顺序为：因此命令流可以输入如下：tb,PIEZ,1tbdata,1,0.171 ! e11tbdata,4,-0.152 ! e21tbdata,7,-0.0187 ! e31tbdata,11,-0.095,0.067 ! e62,e63tbdata,13,0.067 ! e41tbdata,17,0.108,-0.0761 ! e52,e53(5) 网格划分由于本次求解使用的solid226为高阶单元，因此厚度方向划分两层网格即可，即单元长度可选为H的一半0.005m：LESIZE,ALL,5e-3, , , ,1, , ,1Vmesh,1LESIZE命令用来对所选择的线设置网格单元大小，格式为：LESIZE, NL1, SIZE, ANGSIZ, NDIV, SPACE, KFORC, LAYER1, LAYER2, KYNDIV划分后网格如下图：(6) 静态求解输入“/solu”进入求解部分。通过“ANTYPE”命令指定一种分析类型和重启动状态。其格式为：ANTYPE, Antype, Status, LDSTEP, SUBSTEP, Action本例是静态求解，故该命令可为：/soluantype,static !选取静态求解nsel,s,loc,y,0 !定义下表面电压d,all,volt,1000!d,all,ux,0 !约束下表面，若自由边界则不需要nsel,s,loc,y,H1d,all,volt,-1000 !定义上表面电压!d,all,ux,0 !约束上表面，若自由边界则不需要nsel,allsolvefini通过NSEL命令来选择一组节点子集。使用格式为：NSEL, Type, Item, Comp, VMIN, VMAX, VINC, KABS通过D命令在节点上施加DOF约束。其格式为：D, NODE, Lab, VALUE, VALUE2, NEND, NINC, Lab2, Lab3, Lab4, Lab5, Lab6(7) 通用后处理输入“/post1”进入通用后处理部分。本例采用的命令流为：/post1 pldisp,2 !显示变形前的轮廓和变形后的图形plnsol,d,y !显示厚度方向的电流密度save,YshearStaticClamped,dbfiniPLDISP命令用来显示结构的位移，格式为：PLDISP, KUND。PLNSOL命令用来显示节点解。其格式为：PLNSOL, Item, Comp, KUND, Fact, FileID后处理结果如下：l 变形图（位移矢量和）：XZ平面（厚度方向）显然产生了剪切。l 厚度方向电位移图：1.4 结果对比自由边界，厚度方向电位移为：其中，从而夹持边界，厚度方向上的电位移为：由于，所以夹持情况和自由情况下的D2几乎一样大。 从ANSYS计算出的电位移图中可知，即使在几何边界处有类似应力集中的效应，但是整体上石英片厚度方向的电位仪是一致的，大小为C/m2，与理论值十分接近。二、Y切石英板的自由振动2.1 问题描述计算对象为一块石英板，其尺寸为：长度方向（x轴）：L=0.1m厚度方向（y轴）：H=2h=0.01m宽度方向（z轴）：W=0.02m用ANSYS计算石英板在厚度方向振动的模态，并输出前几阶频率。该石英板的材料常数为：密度=2650kg/m3弹性常数矩阵：压电常数矩阵：介电常数矩阵：2.2 理论解控制方程和边界条件为：解出圆频率为：因此频率为：其中，2.3 ANSYS求解(1) 建模与常数输入与Y切石英板静态求解部分一样，这里不再赘述。值得注意的是，由于Y切石英板自由振动为厚度方向的剪切振动，因此厚度方向需要划分更密集的网格。在静态求解时仅划分了两层网格，因此求解厚度方向自由振动时不能再采用静态求解时的网格。本例中对石英板的几何尺寸和网格大小都做出了调整，使得求解所得模态可视性很强，结果也更为准确。(2) 网格划分本例中石英板的几何尺寸为：长度方向（x轴）：L=0.1m厚度方向（y轴）：H=2h=0.01m宽度方向（z轴）：W=0.02m将厚度方向划分为10层，则单元长度为0.001m。由于长度方向和宽度方向比较次要，因此在这两个方向上单元长度可以更大，从而总节点数不至于太大使得求解困难。另外单元也不可过于狭长，因此长度方向和宽度方向的单元长度可以为厚度方向的两倍，即0.002m。APDL如下：lsel,s,length,H1 !选出长度为H1的线段elesize=1e-3LESIZE,ALL,elesizelsel,inve !选出剩下的线段LESIZE,ALL,2*elesizeVmesh,1 !对实体1经行网格划分(3) 模态求解对于模态分析，推荐使用Block Lanczos法（缺省）求解。指令为：/soluantype,modal !选取模态分析modopt,lanb,5 !定义模态分析选项，展开5阶模态mxpand,5(4) 约束与耦合考虑到建模为三维模型，而实际上关注的是Y方向的厚度剪切振动，因此为了使求解效率更高，不产生过多的不需要的模态，这里需要加一些约束与耦合。在厚度方向上的剪切振动中，每一个XZ平面运动都一致，且只能沿着长度方向（X方向）运动，因此约束及耦合如下：nlayer=H1/elesize !厚度方向的层数h=0 !初始平面所在的高度*do,i,1,nlayer !进入循环，i从1到最后一层nsel,s,loc,y,h !选出高度为h的所有节点d,all,uy,0 !约束住y方向上位移d,all,uz,0 !约束住z方向上位移cp,i,ux,all !耦合住x方向上位移h=h+elesize !更新当前平面为下一平面*enddo !结束循环某一层加上约束及耦合后的情况如下：耦合列表如下：耦合共10组，即厚度方向10层全部分别耦合住。(5) 通用后处理后处理我们关心求出的各阶频率及对应的振型，命令流如下：/post1 set,list !列表查看所求模态set,first !查看第一阶模态plnsol,u,sum,2!anmode,10,0.1 !动画查看振型set,next !查看下一阶模态pldisp,2!anmode,10,0.1列表查看求出的频率：第一阶为刚体位移：第二阶开始的振型分别为：N=1阶频率为166302HzN=2阶频率为332245HzN=3阶频率为498518HzN=4阶频率为665262Hz2.4 结果对比(8) 各阶振型图对比ANSYS计算出的各阶振型图如上节所示。理论计算出的振型图如下图所示：二者计算结果一致。(9) 各阶频率对比理论计算各阶频率由下式决定：其中，带入数值得：将ANSYS计算的值与理论解汇入下表：阶数n理论解（Hz）ANSYS解（Hz）误差11660771663020.135%23321543322450.027%34982314985180.058%46643096652620.143%可以看出ANSYS计算出的结果具有十分高的精度。三、（大作业）压电复合梁的静态分析3.1 问题描述如下图所示，一上下对称的压电复合梁。梁厚度方向极化，上下各有一个电极板，板上加有电压。该复合梁长l=0.2m，中间层的弹性梁厚度为he=0.02m，上下两个电极板厚度都为hp=0.005m。本例为一平面应力问题。用ANSYS对其经行静态分析，并与论文”Analysis of Beams with Piezoelectric Actuators”中的结果进行对比。中间层的弹性梁E=2.1E11Pa，=0.3。压电板的物理参数如下：极化方向为z方向极化，但是考虑到ANSYS建二维模型时为XY平面内，因此可将本例中的XZ平面在ANSYS中看做XY平面，即y方向极化。因此各材料常数矩阵如下：弹性常数矩阵为：压电常数矩阵为：介电常数矩阵为：3.2 论文结果论文中计算了两个例子，第一个例子为上下极板电压都为V1=V2=100v，悬臂梁、简支梁对应的中线（z=0）处的扰度曲线为：第二个例子为上极板V1=100v，下极板V2=200v，悬臂梁、简支梁对应的中线处的扰度曲线为：厚度方向应力、应变的变化曲线：中线（z=0）处的位移：3.3 ANSYS求解(1) 建模本例为平面应力问题，可以选择plane13单元。通过选择KEYOPT(1)=7来激活plane13单元的volt自由度。用TB、TBDATA来输入弹性常数矩阵、压电常数矩阵和相对介电常数矩阵时都需要注意ANSYS中的顺序规定。详情可见第一部分，Y切石英的静态求解，这里不再赘述。用于指定单元类型和参数输入的命令流为：!压电单元et,1,plane13,7,2 !2表示平面应力c11=81.3E9c12=0.329E9c13=0.432E9c33=64.5E9c44=25.6E9e31=-5.2e33=15.08e15=0g11=1300g33=1475mp,perx,1,g11 !相对介电常数mp,pery,1,g33tb,ANEL,1 !刚度系数矩阵tbdata,1,c11,c13,c12tbdata,7,c33,c13tbdata,12,c11tbdata,16,c44tb,PIEZ,1 !压电（应力）系数矩阵tbdata,2,e31tbdata,5,e33tbdata,10,e15!弹性单元et,2,plane182mp,EX,2,2.1E11mp,NUXY,2,0.3本例的几何模型十分简单，可以用RECTNG命令直接生成。RECTNG命令的格式为：RECTNG, X1, X2, Y1, Y2因此对应的命令流为：!几何尺寸L=0.2He=0.02Hp=0.005!创建几何体rectng,0,L,-0.5*He-Hp,-0.5*Herectng,0,L,-0.5*He,0.5*Herectng,0,L,0.5*He,0.5*He+Hpaglue,all其中aglue是将生成的三个矩形变成一个整体，从而不会出现多余的面、线、点。(2) 网格划分首先可以用命令ESIZE来制定单元大小。然后要分别对三个矩形进行网格划分。其中下极板对应的面的编号为1，材料编号为1；中间的弹性梁对应的面的编号为4，材料编号为2；上极板对应的面的编号为5，材料编号为1。因此命令流如下：!进行网格划分esize,0.5*Hp !指定单元大小mat,1$type,1amesh,1 !下极板mat,2$type,2amesh,4 !中间层mat,1$type,1amesh,5 !上极板(3) 静态求解静态求解的步骤比较简单。对于本例而言，要注意边界条件的施加。对于悬臂梁，需要约束x=0处的所有结构自由度，命令流如下：!悬臂梁nsel,s,loc,x,0d,all,ux,0d,all,uy,0对于简支梁，需要约束中线（z=0）两端Y方向的位移和左侧端点X方向的位移，命令流如下：!简支梁nsel,s,loc,x,0nsel,r,loc,y,0d,all,ux,0d,all,uy,0nsel,s,loc,x,Lnsel,r,loc,y,0d,all,uy,0对于上下极板的电压的施加类似。最后还要注意极板与弹性梁接触的两个面是接地的，因此电压为0，命令流如下：!定义上极板电压nsel,s,loc,y,0.5*He+Hpd,all,volt,V !V=100!中间层上下面接地nsel,s,loc,y,0.5*Hed,all,volt,0nsel,s,loc,y,-0.5*Hed,all,volt,0!定义下极板电压nsel,s,loc,y,-0.5*He-Hp!注意选择算例case1和case2!d,all,volt,V !case1d,all,volt,2*V !case2(4) 通用后处理简支梁和悬臂梁的位移云图(以case2为例)：为了得出与论文结果相应的几副图，需要将某些路径上的值以图的形式输出。通过命令PATH，PPATH，PDEF和PLPATH可以实现这一点。对于算例1，中线（z=0）上的扰度需要输出为图的形式。其命令流如下： !取出路径z=0path,z0,2ppath,1,0,0ppath,2,L,0!显示该路径上挠度wpdef,w,u,y/axlab,x, x(m) !x轴名称/axlab,y, w(m) !y轴名称plpath,w相应的，悬臂梁和简支梁的图如下：上图：case1，悬臂梁，极值-9.751E-7，论文结果接近但小于-1E-6上图：case1，简支梁，极值2.452E-7，论文结果接近2.4E-7对于算例2，中线（z=0）上的扰度同样需要输出为图的形式：上图：case2，悬