2013有限元单元刚度矩阵(1027讲课用

常用单元的刚度矩阵

r u r r u r =-+= πππεθ22)(2 由于各点在圆周方向上无位移，因而剪应变θr v 和r v θ均为零。将应变写成向量的形式，则{}?? ?? ? ?????? ?????? ???????+??????=??????????????=r w z u z w r u r u rz z r γεεεεθ 根据上式，可推导出几何方程{}[]{})(e B ?ε= 其中几何矩阵[]????????? ?????????? ??= ij ji ki ik jk kj ji ik kj k j i ij kj jk z r z r z r r r r r z r N r z r N r z r N z z z B 000 0),(0),(0),(00021 3.弹性方程和弹性矩阵[D] 依照广义虎克定律，同样可以写出在轴对称中应力和应变之间的弹性方程，其形式为 [])(1 θσσσε+-= z r r u E [])(1 z r u E σσσεθθ+-= [])(1 θσσσε+-=r z z u E rz rz E r τμ)1(2+= 所以弹性方程为{}[]{}εσD = 式中应力矩阵{}{}T rz z r τσσσσθ=

弹性矩阵[]? ? ??????? ???? ?-----+=221000010101)21)(1(μμμμμμμμμμ μμE D 4.单元刚度矩阵[])(e k 与平面问题相同，仍用虚功原理来建立单元刚度矩阵，其积分式为 [][][][]dV B D B k V T e ?=)( 在柱面坐标系中，drdz dV π2= 将drdz dV π2=代入[][][][]dV B D B k V T e ?=)(，则[][][][]rdrdz B D B k T e ??=π2)( 即为轴对称问题求单元刚度矩阵的积分式。与弹性力学平面问题的三角形单元不同，在轴对称问题中，几何矩阵[B]有的元素（如r z r N i ),(等）是坐标r 、z 的函数，不是常量。因此，乘积[][][]B D B T 不能简单地从式 [][][][]rdrdz B D B k T e ??=π2)(的积分号中提出。如果对该乘积逐项求积分，将是一个繁重的工作。一般采用近似的方法：用三角形形心的坐标值代替几何矩阵[B]的r 和z 的值。用[]B 表示在形心),(z r 处计算出的矩阵[B]。其中 3 ) (,3 ) (k j i k j i z z z z r r r r ++= ++= 只要单元尺寸不太大，经过这样处理引起的误差也不大。被积函数又成为常数，可以提出到积分号外面：

ansys质量矩阵刚度矩阵提取

ansys质量矩阵刚度矩阵提取看了这么久了都没人回，查了一些质料终于找到答案了，，下面提供三种方法：方便与其他程序进行接口编程1. Which matrix you would like? element stiffness matrix or full stiffness matrix? element stiffness is within file.emat. full stiffness matrix is within file.full A simple way to dump the matrix is as follow: ------------------- /aux2 fileaux2,file,emat form,long dump,all ------------------- 2. 可以使用/DEBUG命令来得到。详细步骤参见下面的宏文件 finish /clear PI=3.1415926 w1=3 w2=10 w3=6 w4=1.2 r=.8 t=0.08 /PREP7 !* ET,1,SHELL63 R,1,t ET,2,MASS21 R,2,500,500,500,2000,2000,2000,

!* UIMP,1,EX, , ,2e11 UIMP,1,NUXY, , ,0.3, UIMP,1,DAMP, , ,0.2, UIMP,1,DENS, , ,7800, BLC4,0,0,w2,w1 ESIZE,1.5,0, AMESH,all NSEL,S,LOC,X,0.0 D,all, , , , , ,ALL, , , , , allsel,all SFA,all,1,PRES,12 FINISH /OUTPUT,cp,out,, ! 将输出信息送到cp.out文件 /debug,-1,,,1 ! 指定输出单元矩阵 /SOLU SOLVE finish /OUTPUT, TERM ! 将输出信息送到output windows中 ! 这时用编辑器打开cp.out文件，可以看到按单元写出的质量、刚度等矩阵 3. 其原理很简单，即使用ansys的超单元即可解决问题。定义超单元，然后列出超单元的刚度矩阵即可。面是一个小例题，自可明白。 /prep7 k,1 k,2,3000 l,1,2 et,1,beam3 mp,ex,1,2e5 mp,prxy,1,0.3 mp,dens,1,2e3 r,1,5000,2e7,200 lesize,all,,,10 lmesh,all

提取刚度矩阵

============为什么折腾这个文档======== 我有一个计算线性动力学方程组的瞬态、谐响应和静力学的python程序，现希望开发一个将ANSYS组集好的总体矩阵导入该PYTHON程序中的接口。该问题可分解为： [STEP1] [ANSYS]->[包含矩阵信息的文件] [STEP2] [包含矩阵信息的文件]->[python通用数据对象] [STEP3] [python通用数据对象]->[程序特定数据对象]->[进行计算] 因此检索了一些帖子，基本上完成了这项工作，本文是对[STEP1]和[STEP2]的整理，并且利用[STEP3]对结果进行了验证 ============主要内容================== 1，了解从ANSYS中提取总体矩阵和载荷向量的方法； 2，了解提取出来的矩阵是怎样表示的； 3，说明在Python中，如何读取这样的矩阵； 4，构造一个简单的算例，说明整个【建模】-【提取】-【读取】过程及其正确性； =========站内检索综述==================== 检索词：提取矩阵得到21个结果，代表性的帖子有下面这9个：编号[1] 标题：ansys中怎样提取质量，刚度，阻尼矩阵？地址：https://www.360docs.net/doc/9e14005361.html,/forum-vi ... fromuid-159019.html 要点：pengweicai给出了一段网上最常见的提取代码,该程序以fortran 写成，可以利用.full文件以及一些列约定将ANSYS中的总体矩阵读入FORTRAN中。编号[2] 标题：如何得知HBMAT命令提取的质量、刚度矩阵对应的自由度？地址：https://www.360docs.net/doc/9e14005361.html,/forum-vi ... fromuid-159019.html 要点：提出了使用HBMAT命令提取稀疏矩阵时常见的问题：我们如何知道提取出来的信息是怎么储存的呢？编号[3] 标题：[分享]ANSYS中整体、单元刚度和质量矩阵的提取地址：https://www.360docs.net/doc/9e14005361.html,/forum-vi ... fromuid-159019.html 要点：在该帖子的7楼，其实已经给出了帖子[2]中问题的解答，即HBMAT 中提取出来的矩阵是Harwell-Boeing格式的，并且给出了该格式的细节，可惜是英文的，没引起多少关注。编号[4] 标题：帮我看看提取的刚度与质量矩阵地址：https://www.360docs.net/doc/9e14005361.html,/forum-vi ... fromuid-159019.html 要点：这个帖子所示的矩阵并非是使用HBMAT命令提出出来的，而应该是SELIST命令列举出来的未压缩的矩阵，后续楼层的回帖给了大家一个提示，即有可能提取出来的矩阵是引入了边界条件的（即删除了被约束的行和列的）。编号[5] 标题：提取刚度矩阵的问题地址：https://www.360docs.net/doc/9e14005361.html,/forum-vi ... fromuid-159019.html 要点：本帖作者的工作是基于单元刚度矩阵的，因此ANSYS中提取的单元刚度矩阵是否处于总体坐标系就成为问题。该问题并非本文内容，但仍值得关注。编号[6] 标题：提取刚度矩阵丢失节点的问题地址：https://www.360docs.net/doc/9e14005361.html,/forum-vi ... fromuid-159019.html 要点：帖子[5]作者的又一帖，在这里帖子[5]的问题得到了欧阳中华老师的回答。编号[7] 标题：提取刚度矩阵的ANSYS操作过程地址：https://www.360docs.net/doc/9e14005361.html,/forum-vi ... fromuid-159019.html 要点：实际上这就是使用HBMAT从ANSYS中提取总体矩阵的全过程！只是还有一些细节待确定。编号[8] 标题：提取整体刚度矩阵、质量矩阵及阻尼矩阵的简单方法地址：https://www.360docs.net/doc/9e14005361.html,/forum-vi ... fromuid-159019.html 要点：给出了利用“不减缩的”子结构方法来得到总体矩阵的方法（这也是网络上常见的代码之一）编号[9] 标题：质量矩阵、刚度矩阵如何提取？地址：https://www.360docs.net/doc/9e14005361.html,/forum-vi ... fromuid-159019.html 要点：16443在5楼的回帖中给出了提取刚度矩阵的三种方法 =======站外检索略述======================== 百度检索：提取矩阵比较好的帖子有：编号[10] 来源:百度文库标题:怎样从ansys中提取单元刚度矩阵与质量矩阵地址:https://www.360docs.net/doc/9e14005361.html,/view/3cf5e567f5335a8102d220d9.html 要点:这应该就是16443在帖子[9]中回复的内容了，全面的总结了在帖子[3,4,5,9]中涉及的问题。编号[11] 来源:中华钢结构标题:ansys刚度矩阵Harwell-Boeing格式的具体含义讨论地址:https://www.360docs.net/doc/9e14005361.html,/forum/viewthread.php?tid=184007 要点:如题，后续楼层给出了一些将矩阵读入ANSYS的APDL（好不容易读出来，又读进去干嘛呢……）编号[12] 来源:simwe 标题:关于ANSYS（质量、刚度、阻尼）矩阵Harwell-boeing格式数据的说明地址:https://www.360docs.net/doc/9e14005361.html,/archiver/tid-924778.html 要点:比[11]更透彻的HB格式说明! ============================================================= =======1.从ANSYS中提取总体矩阵的方法================================= ============================================================= 1,用/DEBUG命令 2,子结构法

一般单元在局部坐标系下的单元刚度矩阵

9.3 一般单元在局部坐标系下的单元刚度矩阵 1.杆端内力与位移关系回顾（轴向）；；（弯曲）； 2.公式推导（图1）图1 杆件性质：长度l，截面面积A，截面惯性矩I，弹性模量E；杆端位移u、v、θ。（1）（2）列成矩阵形式：

（3）即：（4）局部坐标系下单元刚度矩阵：（5） 9.4 梁单元 1.简支梁简支梁单元见图1。图1 说明：（a）梁单元通常忽略轴向变形；（b）图10-3中；相应的力分量也应该为零；（c）依据刚度矩阵的物理意义，可以由一般单元的刚度矩阵生成梁单元矩阵。即去掉位移分量为零的相应行和列。

即：单元刚度方程：单元刚度矩阵：（1） 2.悬臂梁等思考：建立图2的单元刚度矩阵：（固定端位移为零；自由端有转角和竖向位移) 图2 图a：图b： 3.桁架仅有轴向位移 9.5 单元刚度系数的物理意义 1.单元刚度系数的意义一般地，第j 个杆端位移分量取单位值1，其它杆端位移为0 时所引起的第i个杆端力分量的值。

例：的物理意义：当第3个杆端位移分量时引起的第5个杆端力分量。对称性（反力互等定理） 3.奇异性（，不存在逆矩阵）根据式可由杆端位移求解杆端力，且是唯一解。但由杆端力求杆端位移，可能无解，如有解也是非唯一解。说明：已知6个杆端力分量，（a）无法保证力状态的合法性——可能造成无解；（b）无法确定杆的支承条件——可能造成非唯一解。 9.6 单元坐标转换矩阵的物理意义 1.问题的提出单元刚度矩阵——单根杆；多根根组成的复杂结构呢？（图1）

图1 分析（a）从数学的角度理解整体坐标系（xy）与局部坐标系（）的区别；（b）力分量应向整体坐标系转换，图f给出了两种坐标系下力分量之间的数学关系：。同理： 2.公式推导矩阵形式：（1）同理：（2）

7.4 单元刚度矩阵组装及整体分析

根据全结构的平衡方程可知，总体刚度矩阵是由单元刚度矩阵集合而成的个结构的计算模型分成个单元，那么总体刚度矩阵可由各个单元的刚度矩阵组装而成，即是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的将总体坐标轴分别用表示，对某单元有式中，和分别是局部坐标系和总体坐标系下的单元结点位移向量是该单元在总体坐标系下的单元刚度矩阵

. 将单元结点的局部编号换成总体编号，

其中右上角的上标表示第单元所累加上的子矩阵具有相同的下标，的那些子矩阵的累加总体刚度矩阵第行的非零子矩阵是由与结点相联系的那，从环绕点各单元移置而来的结点载荷为式中表示对环绕结点的所有单元求和，环绕结点的各单元施加于结点的结点力为

.因此，结点的平衡方程可表示为得到以结点位移表示的结点的平衡方程，为整体刚度矩阵，为全部结点位移组成的向量，为全部结点载荷组成的向量式中，是总体坐标系下的结点载荷向量，为坐标转换阵 . 构是处于自由状态，在结点载荷的作用下，结构可以产生任意的刚体位移的条件下，仍不能通过平衡方程惟一地解出结点位移.

约束的种类包括使某些自由度上位移为零，，或给定其位移值，还有给定支承刚为了理解这个方法，我们把方程分块如下：其中，假设是给定的结点位移；是无约束的（自由）结点位移因而是已知的结点力；

其中，不是奇异的，因而可以解方程（一旦知道了，求得未知结点力. 殊情况下，我们可以删除对应于的各行和各列（即删行删列法），故可把方程简写为由于全部给定的结点位移通常都不能在位移向量的开始或终了，故分块法的编号方法是很麻烦因此，为了引入给定的边界条件，可以采用下述等价的方法如果把给定为，则载荷向量为结点自由度总数

提取单元刚度矩阵

单元刚度矩阵的提取刚度矩阵在有限元求解过程中扮演者非常重要的角色，以最小位能原理求解过程为例最终越是转换为含有结构刚度矩阵的能量泛函的取值问题。有限元过程中涉及到三类刚度：单元刚度矩阵，组合结构刚度矩阵和最终求解刚度矩阵。其中单元刚度矩阵：仅与单元的自身自由度有关，同一编号的单元矩阵的维数是固定。组合结构刚度：矩阵根据求解的初始变量个数决定刚度矩阵的维数，属于单元组装后的初始刚度，维数和整个单元初始变量个数相等。最终求解刚度矩阵：代入边界条件简化后的刚度。以《Finite Element Analysis-Theory and Application With ANSYS》中的梁单元例子为例，解释刚度提取过程：此模型的单元刚度矩阵：（学则beam3梁单元后，该单元包含两节点，每个节点具有三个自由度，因此对应单元刚度矩阵为6*6的方阵）

组合结构刚度矩阵：（该结构含有三个节点，每个节点具有三个原始自由度，因此组合结构刚度矩阵具有9*9阶的形式）最终求解刚度矩阵：（由于边界条件的存在，该结构中，1,3点的自由度不存在，求解参数中有六个参数已知，因此对最终求解刚度矩阵为三阶方阵）通过最终的刚度矩阵组成的方程，求解出2节点的位移解，再以这些原始解得出应力，应变，支反力的其他的解。 ansys实现过程：提取思路如下：通过/debug提取单元刚度矩阵，通过filname.full文件提取后两者的矩阵 ansys实现过程如下： finish /clear /filname,k,1 /prep7 N,1 N,2,120 N,3,120,-108 et,1,beam3 mp,ex,1,3.0e7 mp,prxy,1,0.3 R,1,7.65,204,10 E,1,2 E,2,3 /debug,-1,,,1,,,,,

有限元单元刚度矩阵

实现功能：输入欲求的单元号,得到该单元的面积和该单元刚度矩阵；输入完所有的单元号,得到所有的单元刚度矩阵后,直接得到半带宽存储的数组。（所有的单元刚度矩阵和半带宽数组分别存在所有单元的刚度矩阵.txt和SK矩阵.txt）第1单元面积为：0.5000 弹性模量、泊松比和厚度分别为:100.0000 0.3000 0.1000 单元1单元的应力矩阵 -109.8901 -32.9670 109.8901 0.0000 0.0000 32.9670 -32.9670 -109.8901 32.9670 0.0000 0.0000 109.8901 -38.4615 -38.4615 0.0000 38.4615 38.4615 0.0000 第1单元的单元刚度阵为 7.4176 3.5714 -5.4945 -1.9231 -1.9231 -1.6484 3.5714 7.4176 -1.6484 -1.9231 -1.9231 -5.4945 -5.4945 -1.6484 5.4945 0.0000 0.0000 1.6484 -1.9231 -1.9231 0.0000 1.9231 1.9231 0.0000 -1.9231 -1.9231 0.0000 1.9231 1.9231 0.0000 -1.6484 -5.4945 1.6484 0.0000 0.0000 5.4945 半带宽存储的数组 7.4176 3.5714 -1.9231 -1.6484 -5.4945 -1.9231 7.4176 -1.9231 -5.4945 -1.6484 -1.9231 0.0000 9.8901 -2.3810 -4.9451 4.7619 -3.0220 -0.4579 9.8901 4.7619 -4.9451 -0.7326 0.5495 0.0000 20.8791 -2.3810 0.5495 -2.3810 -10.989 1.6484 13.7363 -2.3810 -3.0220 1.9231 -3.8462 0.0000 3.4341 1.1905 -0.9615 1.9231 0.0000 0.0000 5.2198 1.6484 -2.7473 0.0000 0.0000 0.0000 11.9506 -3.5714 0.0000 0.0000 0.0000 0.0000 6.5934 0.0000 0.0000 0.0000 0.0000 0.0000 #include #define jds 5 #define d 2 main() {int nj,ne,i,j,k,a,n,*lnd; int JDS[3];xj,yj,xo,yo; float B[3][6],b1[3][6];D[3][3],E,u,t;S[3][6],s1[6][3];K[6][6];Kk[2*jds][2*jds]={{0}}; float SK[2*jds][(d+1)*2]={{0}}; float x1,x2,x3,y1,y2,y3,ae,*xy; FILE *fp1,*fp2,*fp3,*fp4,*fp5,*fp7; fp1=fopen("节点坐标.txt","r"); fp2=fopen("单元节点编号.txt","r");fp3=fopen("材料参数.txt","r"); fp4=fopen("输出结果.txt","a"); fp5=fopen("节点总数和单元总数.txt","r");fp7=fopen("SK矩阵.txt","w"); fscanf(fp5,"%d",&nj); fscanf(fp5,"%d",&ne); xy=(float*)malloc(nj*2*sizeof(float)); lnd=(int*)malloc(ne*3*sizeof(int));

MatLab在有限元刚度矩阵推导中的应用

第２６卷第２期重庆交通学院学报ｖ０１．２６Ｎｏ．２２００７年４月ＪＯＵＲＮＡＬＯＦＣＨＯＮＧＱＩＮＧ：姒０７ｒＯＮＧＵＮⅣＥＲＳｒｒＹＡｐｒ．．２００７ ?＿＿＿＿Ｉ＿＿＿＿目一ｒ－－－！ｒ自目ｌ＿＝＝日日一ＭａｔＬａｂ在有限元刚度矩阵推导中的应用周水兴１’２，陈山林１（１．重庆大学，重庆４０００４５；２．重庆交通大学土木建筑学院，重庆４０００７４）摘要：ＭａｔＬａｂ具有强大的符号运算功能．以平面梁单元刚度矩阵为例，详细介绍ＭａｔＬａｂ在单元刚度矩阵推导中的具体应用，编写了平面梁单元刚度符号运算程序，运算结果与手工推导完全一致，该方法可进一步推广到其它单元甚至到更复杂非线性单元刚度矩阵推导中．关键词：ＭａｔＬａｂ；符号运算；刚度矩阵；推导中图分类号：０２４２．２１文献标识码：Ａ文章编号：１００１．７１６Ｘ（２００７）０２．００２９．０３ＭａｔＬａｂＡｐｐｌｉｃａｔｉｏｎｔｏｔｈｅＭａｔｒｉｘＤｅｄｕｃｔｉｏｎｏｆＦｉｎｉｔｅＥｌｅｍｅｎｔＳｔｉｆｆｎｅｓｓＺＨＯＵＳｈｕｉ—ｘｉｎ９１．．。ＣＨＥＮＳｈａｎ．１ｉｎｌ（１．ＣｈｏｎｇｑｉｎｇＵｎｉｖｅｒｓｉｔｙ，Ｃｈｏｎｇｑｉｎｇ４０００４５，Ｃｈｉｎａ；２．ＳｃｈｏｏｌｏｆＣｉｖｉｌＥｎｇｉｎｅｅｒｉｎｇ＆Ａｒｃｈｉｔｅｃｔｕｒｅ，ＣｈｏｎｇｑｉｎｇＪｉａｏｔｏｎｇＵｎｉｖｅｒｓｉｔｙ，Ｃｈｏｎｇｑｉｎｇ４０００７４，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：ＭａｔＬａｂｈａｓａｓｔｒｏｎｇｓｙｍｂｏｌｉｃｏｐｅｒａｔｉｏｎｆｕｎｃｔｉｏｎ．Ｗｉｔｈｔｈｅｅｘａｍｐｌｅｏｆｓｔｉｆｆｎｅｓｓｍａｔｒｉｘｏｆｐｌａｎｅｂｅａｍｅｌｅｍｅｎｔ，ｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｄｅｄｕｃｉｎｇｅｌｅｍｅｎｔｓｔｉｆｆｎｅｓｓｍａｔｒｉｘｉｓｉｎｔｒｏｄｕｃｅｄｉｎｄｅｔａｉｌｆｏｒｔｈｅｆｉｒｓｔｔｉｍｅ，ａｎｄｔｈｅｓｙｍｂｏｌｉｃｏｐｅｒａｔｉｏｎｐｒｏｃｅｄｕｒｅｏｆｐｌａｎｅｂｅａｍｅｌｅｍｅｎｔｉｓｄｅｖｄｏｐｅｄ，ａｎｄｔｈｅｏｐｅｒａｔｉｏｎｒｅｓｕｌｔｉｓｃｏｍｐｌｅｔｅｌｙｉｎａｃｃｏｒｄａｎｃｅｗｉｔｈｔｈｅｈａｎｄｉｃｒａｆｔｄｅｄｕｃｔｉｏｎ．Ｔｈｉｓｍｅｔｈｏｄｃａｎｆｕｒｔｈｅｒｅｘｐａｎｄｔｈｅａｐｐｌｉｃａｔｉｏｎｏｆｄｅｄｕｃｉｎｇｍａｔｒｉｘｔｏｏｔｈｅｒｅｌｅｍｅｎｔｓｅｖｅｎｍｏｒｅｃｏｍｐｌｉｃａｔｅｄｎｏｎｌｉｎｅａｒｅｌｅｍｅｎｔ．ＫｅｙｗｏｒｄｓｉＭａｔｌ．ａｂ，ｓｙｍｂｏｌｉｃｏｐｅｒａｔｉｏｎ，ｅｌｅｍｅｎｔｓｔｉｆｆｎｅｓｓｍａｔｒｉｘ，ｄｅｄｕｃｅＭａｔＬａｂ是美国ＭａｔｈＷｏｒｋ公司开发的用于数值计算、概念设计、算法研究、建模仿真、实时实现理想的集成环境，因其完整的专业体系和强大的运算功能，已广泛应用于工业、电子、信号处理、控制、医疗、建筑、教学等各个领域…．ＭａｔＬａｂ本意是矩阵实验室（ＭａｔｒｉｘＬａｂｏｒａｔｏｒｙ），以向量和矩阵运算为基础，强大的数值计算功能是其最具代表性的特点．ＭａｔＬａｂ的数学计算包括数值计算和符号计算，前者主要用于矩阵、向量的各种数值计算；后者主要是计算式中带有符号变量、表达式的运算．该运算符由符号数学工具箱提供支持，有复合、简化、微分、积分及求解代数方程和微分方程的工具．ＭａｔＬａｂ符号数学工具箱中的工具建立在功能强大的Ｍａｐｌｅ软件基础上．有限元是近代数值计算最有效方法之一．有限元法的基础是单元刚度矩阵的推导，目前，刚度矩阵推导已有一个相对固定的模式，而烦琐、复杂的矩阵运算、微分、积分是刚度矩阵推导中的主要内容．通常，这种矩阵运算是由手工来完成的，工作量大，而且极易出错．利用ＭａｔＬａｂ丰富的符号运算功能，完成刚度矩阵推导过程中的符号运算，不但速度快，而且准确性高．Ｍａｆｌａｂ的命令和函数相当多，本文仅就涉及到ＭａｔＬａｂ符号运算作一些简要介绍，然后结合平面梁单元刚度矩阵的推导，用ＭａｔＬａｂ完成该单元刚度矩阵的符号运算，结果完全一致．１ＭａｔＬａｂ符号运算的一般规定１．１符号标变量和表达式的创建收稿日期：２００６－０２－０７；修订日期：２００６－０２－２３基金项目：重庆市教委项目（ＫＪ０５０４０７）；交通部西部交通建设科技项目（２００３３１８８１４２０２４）作者简介：周水兴（１９６７－），男，浙江嘉兴人，教授，主要从事桥梁设计理论研究．ｅ－ｍａｌｌ：ｓｈｕｉｘｉｎｇｚｈｏｕ＠ｃｑｕｃ．ｅｄｕ．ｃｎ

基于matlab的有限元法分析平面应力应变问题刘刚

姓名：刘刚学号：15 平面应力应变分析有限元法 Abstruct:本文通过对平面应力/应变问题的简要理论阐述，使读者对要分析的问题有大致的印象，然后结合两个实例，通过MATLAB软件的计算，将有限元分析平面应力/应变问题的过程形象的展示给读者，让人一目了然，快速了解有限元解决这类问题的方法和步骤！一.基本理论有限元法的基本思路和基本原则以结构力学中的位移法为基础，把复杂的结构或连续体看成有限个单元的组合，各单元彼此在节点出连接而组成整体。把连续体分成有限个单元和节点，称为离散化。先对单元进行特性分析，然后根据节点处的平衡和协调条件建立方程，综合后做整体分析。这样一分一合，先离散再综合的过程，就是把复杂结构或连续体的计算问题转化简单单元分析与综合问题。因此，一般的有限揭发包括三个主要步骤：离散化单元分析整体分析。二.用到的函数 1. LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym,p) (K k I f) (k u) (k u A) (E NU t) 三.实例例1.考虑如图所示的受均布载荷作用的薄平板结构。将平板离散化成两个线性三角元，假定E=200GPa，v=，t=0.025m,w=3000kN/m. 1.离散化 2.写出单元刚度矩阵

通过matlab 的LinearTriangleElementStiffness 函数，得到两个单元刚度矩阵1k 和2k ，每个矩阵都是6 6的。 >> E=210e6 E = >> k1=LinearTriangleElementStiffness(E,NU,t,0,0,,,0,,1) k1 = +006 * Columns 1 through 5 0 0 0 0 0 0 0 0 Column 6 >> NU= NU = >> t= t = >> k2=LinearTriangleElementStiffness(E,NU,t,0,0,,0,,,1)

ansys质量矩阵刚度矩阵提取说课材料

a n s y s质量矩阵刚度矩阵提取

第12章模态分析21

第12章模态分析 12.1 模态分析概述模态分析是ANSYS 中分析结构自然频率和模态形状的方法；它假设：①结构刚度矩阵和质量矩阵不发生改变；②除非指定使用阻尼特征求解方法，否则不考虑阻尼效应；③结构中没有随时间变化的载荷。在无阻尼系统中，结构振动方程如下 []{}[]{}{}0=+u K u M （12-1）式中，[]M 为质量矩阵；[]K 为刚度矩阵；{}u 为节点加速度向量；{}u 节点位移向量。其中刚度矩阵可以包括预应力效应带来的附加刚度。对线性系统而言，自由振动满足下面方程 {}{}t u i i ω?cos = （12-2）式中，{}i ?为第i 阶模态形状的特征向量； i ω第i 阶自然振动频率；t 时间。将（12-2）代入方程（12-1），得到 [][]() {}{}02=+-i i K M ?ω （12-3）从式（12-3）中得到结构的振动特征方程为 [][]02=+-K M i ω （12-4）通过式（12-4）可以求出第i 阶自然振动频率i ω，进而代入（12-3）可以求出第i 阶模态形状的特征向量{}i ?。将{}i ?对质量矩阵[]M 进行归一化处理，使用命令MODOPT,,,,,,,OFF ，可以得到 {}[]{}1=i T i M ?? （12-5）如果{}i ?，向自身做归一化处理，使用命令MODOPT ,,,,,,ON ，那么{}i ?中最大的向量坐标将归一化为1.0。如果使用缩减模态提取方法，使用MODOPT,REDUC ，第i 阶模态形状的特征向量{}i ?可以通过使用MXPAND 命令进行扩展。

12.2模态分析过程 ANSYS的模态分析是线性分析的一种，对于任何非线性特性，如塑性和接触（间隙）单元，在模态分析中将被忽略。模态分析过程由4个主要步骤组成，即前处理、加载与求解、扩展模态，以及查看结果和后处理。 12.2.1前处理建模是指建立分析的有限元数学模型，包括建立几何模型和划分网格，模态分析的建模过程与一般的建模过程并没有实质性的区别，具体建模可以参见第三章。但根据模态分析的特点，需要注意以下几点： ?定义材料特性时，必须考虑质量的问题。如果最终得到的模型中没有任何质量，那么质量矩阵将为[]0，而无法求解系统的固有频率。 ?模态分析只考虑材料的线性行为。材料可以为线性各向同性、正交各向异性、温度无关和温度有关等类型，必须定义材料的杨氏模量和质量相关属性。对于可能定义的非线性特性，ANSYS在求解时都将忽略。模态分析只考虑网格单元的线性行为，对于非线性的单元类型将会被视为线性单元处理，例如在结构中定义了接触单元，在分析中将计算接触单元初始状态的刚度矩阵，而将此刚度矩阵应用到分析的其他任何时候。对于预应力分析，模态分析将接触单元的刚度矩阵取为静态预应力分析结束时的刚度矩阵。如果定义特殊的阻尼单元类型（如COMBIN14, COMBIN37等），必须按单元的要求定义需要的实常数。 12.2.2加载与求解在这个步骤中要定义分析类型和分析选项，施加载荷，指定加载阶段选项，并进行固频率的有限元求解。应在求解前设置模态扩展选项，或在得到初始解后，对模态进行扩展以供查看。1．设置分析类型首先进入求解器，并使用ANTYPE命令或GUI交互的方式，定义求解类型为模态分析。具体操作方法如下。命令方式： ANTXPE, 2 GUI方式：选择Main Menu > Solution > Analysis Type > New Analysis命令，弹出New Analysis对话框，在对话框中选中Modal，单击OK按钮确认。 2．设置分析选项

第三节刚度矩阵(汇编)

第三节刚度矩阵 ——节点载荷与节点位移之间的关系一、单元刚度矩阵 1. 单元刚度矩阵 xj 单元e 是在节点力作用下处于平衡。节点i 的节点力为 {}T i xi yi R R R ??=?? （i , j , m 轮换）则单元e 的节点力列阵为 {} T e T T T m i j T xm ym xi yi xj yj R R R R R R R R R R ??? ? ???? = = 单元应力列阵为 {} T e x y xy σσστ???? =

假定弹性体的所有节点都产生一虚位移，单元e 的三个节点的虚位移为 {} * ***** *e T m m i i j j u v u v u v δ??? ? = 单元虚应变列阵为 {} ****T x y xy εεεγ???? ?? = 参照式（3-7），则单元虚应变为 {} {}**e e B εδ????= 作用在弹性体上的外力在虚位移上所做的功为： {}{}* e T e R δ?? ?? ? 单元内的应力在虚应变上所做的功为： {}{}* T e tdxdy εσ? ?? ?? ? ?? 根据虚位移原理，可得单元的虚功方程 {} {}{} {}**e T T e e R tdxdy δεσ? ???? ? ??? ?? = ?? 或 {}{} {}{}* * e T T T e e B R tdxdy δδσ? ? ????? ? ??? ? ? ? ? =??

故有 {} {}e T B R tdxdy σ? ???? = ?? 将式（3-10）代入，的 {} {}{}e e e T T D B D B R B B tdxdy tdxdy δδ?? ??? ???????????? ?????????== ???? （3-27）简记为 {}{}e e e k R δ?? ?? = （3-29） --------上式表征单元节点力与节点位移之间的关系，称为单元刚度方程（单元平衡方程）其中 T e D B B k tdxdy ? ????? ??????????? = ?? （3-28） e k ????称之为单元刚度矩阵（简称为单刚），是66?矩阵。如果单元的材料是均质的，矩阵D ????中的元素也是常量，且在三角形常应变的情况下，矩阵B ????中的元素也是常

单元刚度矩阵(等参元)MATLAB编程

《有限元法》实验报告专业班级力学（实验）1601 姓名田诗豪学号 10 提交日期

实验一（30分）一、实验内容编写一个计算平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的MATLAB 函数文件[B3,S3,K3] = ele_mat_tri3(xy3,mat)，其中：输入变量xy3为结点坐标数组，mat为材料参数矩阵；输出变量B3为应变矩阵，S3为应力矩阵，K3为单元刚度矩阵。（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码通用函数 function [B3,S3,K3] = ele_mat_tri3(xy3,mat) %生成平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的功能函数 %*********变量说明**************** %xy3------------------结点坐标数组 %mat------------------材料参数矩阵（弹性模量，泊松比，壁厚） %B3-------------------应变矩阵 %S3-------------------应力矩阵 %K3-------------------单元刚度矩阵 %********************************* xyh=[1,xy3(1,1),xy3(1,2);1,xy3(2,1),xy3(2,2);1,xy3(3,1),xy3(3,2)]; A=*det(xyh); A=abs(A); D=mat(1)/(1-mat(2)^2)*[1,mat(2),0;mat(2),1,0;0,0,(1-mat(2))/2]; b=zeros(1,3);c=zeros(1,3); %********************************* for i=1:3 if i==1 j=2; m=3; elseif i==2 j=3; m=1; else j=1; m=2; end b(i)=xy3(j,2)-xy3(m,2); c(i)=xy3(m,1)-xy3(j,1);

ansys 常见技巧汇总

一、前处理 1.实体显示*.sat、*.x_t等外部导入模型/facet,fine /replot Gui: Utility Menu>PlotCtrls>Style>Solid Model Facets 2.修改ansys背景jpgprf，500，100，1 /replot 3.隐藏坐标系的显示/triad,off /replot Gui: Utility Menu>PlotCtrls>Window Controls>Reset Window Options Utility Menu>PlotCtrls>Window Controls>Window Options 4.设置参考温度TREF, TREF Gui：Main Menu>Solution>Define Loads>Settings>Reference Temp 5.显示单元实际形状/eshape,1.0 Gui: Utility Menu>PlotCtrls>Style>Size and Shape 6.透明显示单元、体、面/TRLCY, Lab, TLEVEL, N1, N2, NINC Gui: Utility Menu>PlotCtrls>Style>Translucency 7.显示编号/PNUM, Label, KEY Gui: Utility Menu>PlotCtrls>Numbering 8.导入hypermesh有限元模型/input,filename,prp Gui: Utility Menu>File>Read Input from 9.导入abaqus格式的有限元模型/input,filename,inp Gui：Gui: Utility Menu>File>Read Input from 10.ansys作为fluent前处理输出 cdwrite,db,filename,cdb gui: Main Menu>Preprocessor>Archive Model>Write 11.不显示单元轮廓线 /gline,1,-1 Gui: Utility Menu>PlotCtrls>Style>Edge Options 12.显示施加到几何元素上的约束 dtran /replot Gui：Main Menu>Preprocessor>Loads>Define Loads>Operate>Transfer to FE>Constraints 13.显示施加到几何元素上的面载荷 sftran /replot Gui:Main Menu>Preprocessor>Loads>Define Loads>Operate>Transfer to FE>Surface Loads 14.显示载荷标记及数值 /pbc,f,,2 Gui: Utility Menu>PlotCtrls>Symbols 15.设置显示容差BTOL, PTOL默认值PTOL为1e-5，可以根据需要修改GUI: Main Menu>Preprocessor>Modeling>Operate>Booleans>Settings 16.如何使用用户定义用户自定义矩阵 Matrix 27用户定义用户自定义矩阵，由单元选项控制定义质量、刚度或阻尼矩阵，你只要在同一组接点，分别定义三次MATRIX27单元（KEYOPT（2）分别为2，4，5）即可，然后在定义实常数时，分别定义三种单元对应的质量、刚度、阻尼矩阵系数。 17.ANSYS的UNDO功能 ANSYS的UNDO功能，多数人都认为ansys没有undo功能。其实这个功能一直就存在，在安装目录\apdl\start100.ans(10.0版，其他版本相应数值变化)，后面加上两行命令 /undo,on$*abbr,undo,undo.启动ansys以后就会出现一个undo的快捷工具18.运算完成后电脑自动关机喜欢用apdl的朋友可能会碰到这么一个麻烦：就是当运算量较大的时候不知道