ANSYS的仿真及教学日记第章Current

1、ANSYS的仿真及教学日记1.有限元的基础知识21.1.函数逼近的思想21.2.差分法求解偏微分方程41.3.差分法求解偏微分方程51.4.等参变换、变分法、单元刚度矩阵相关61.5.MATLAB有限元编程求解力学问题91.6.MATLAB直接使用FEM求解Excel中的简单温度场分布151.7.有限元的线性方程组的求解162.有限元建模182.1.ANSYS10.0的起始菜单的说明：182.2.ANSYS仿真中常用的单元：193.静力学分析253.1.壳体结构的有限元分析253.2.几个教学案例的分析要点：253.3.结构分析线性方程求解的基本算法291. 有限元的基础知识1.1. 函数逼近

2、的思想(a) 基于全域x0,xL的函数逼近（全域高阶插值或拟合的思想）(b) 基于子域xi,xi+1的函数逼近（子域低阶分段插值或分段拟合的思想-有限元的思想）1.2. 差分法求解偏微分方程Excel演示：上边（100,左为75，右为50，底为0）的温度场1.3. 等参变换、变分法、单元刚度矩阵相关#.直接刚度矩阵建立与求解1) 什么是应力、应变、位移2）什么是插值函数？什么是高次单元？什么是低次单元？3）直接刚度法？4）有限元法求解的是什么方程？#.刚度与刚度矩阵弹簧的变形与力的关系k是弹簧系数，是弹簧的固有特性，和截面积、长度、形状有关系变成矩阵的形式：也就是对于杆件：，节点位移，节点外力

3、该节点的内力为：（正负号是考虑到方向）而内力和外力的平衡有，, 因此该方程就是刚度矩阵#.等参变换由于子单元的形状复杂多变，采用等参变换后的坐标的刚度矩阵除了一些非常简单的情况，基本上不能进行显式积分，而需要进行数值积分，如采用Newton-Cotes积分等等，基于对积分方程的各种假设，比如使用二次的线性来代替，或者三次线性函数代替，或者其他假设的函数等等，这样的单元刚度矩阵有各种各样的方法了。3. 单元的各个参量的插值函数这是一个自由度的情况，该单元位移场为:取首两项作低阶插值：节点条件为：带入得到，带入得到：N叫做形状函数矩阵，具体为q是节点位移矩阵#. 整体坐标系和局部坐标系#. 能量原

4、理（虚位移原理、最小势能原理）、变分法和有限元的建立物理方程就是微分方程加上边界条件，对微分方程的原方程的求解是非常困难，甚至不可能的，变分法就是先假设一个泛函（这个泛函是容易处理的标准函数，如线性函数等），泛函的求导就是这个微分方程，通过泛函的求极值来满足边界条件，最后还是得到单元刚度矩阵。通过胡克定律和能量方程等可以求出中性层的挠度方程为 （p为压强）最小势能原理得到（假设一个正弦函数为基函数）：应变能为：W=1.4. MATLAB有限元编程求解力学问题四杆桁架结构的有限元分析(Bar2D2Node) 如图所示的结构，各个杆的弹性模量和横截面积都为, 。试基于MATLAB平台求解该结构的节

5、点位移、单元应力以及支反力。图 四杆桁架结构解答：对该问题进行有限元分析的过程如下。 （1） 结构的离散化与编号对该结构进行自然离散，节点编号和单元编号如图。（2）计算各单元的刚度矩阵(基于国际标准单位) 建立一个工作目录，将所编制的用于平面桁架单元分析的4个MATLAB函数放置于该工作目录中，分别以各自函数的名称给出文件名，即：Bar2D2Node_Stiffness，Bar2D2Node_Assembly，Bar2D2Node_Stress，Bar2D2Node_Forces。然后启动MATLAB，将工作目录设置到已建立的目录中，在MATLAB环境中，输入弹性模量E、横截面积A，各点坐标x

6、1,y1,x2,y2,x3,y3,x4,y4,角度alpha 1, alpha 2和alpha 3,然后分别针对单元1，2，3和4，调用4次Bar2D2Node_Stiffness，就可以得到单元的刚度矩阵。相关的计算流程如下。E=2.95e11;A=0.0001;x1=0;y1=0;x2=0.4;y2=0;x3=0.4;y3=0.3;x4=0;y4=0.3;alpha1=0;alpha2=90;alpha3=atan(0.75)*180/pi;k1=Bar2D2Node_Stiffness (E,A,x1,y1,x2,y2,alpha1) k1 = 73750000 0 -73750000

7、0 0 0 0 0 -73750000 0 73750000 0 0 0 0 0k2=Bar2D2Node_Stiffness (E,A,x2,y2,x3,y3,alpha2) k2 = 1.0e+007 * 0.0000 0.0000 -0.0000 -0.0000 0.0000 9.8333 -0.0000 -9.8333 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -9.8333 0.0000 9.8333k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3,alpha3) k3 = 1.0e+007 * 3.7760 2.83

8、20 -3.7760 -2.8320 2.8320 2.1240 -2.8320 -2.1240 -3.7760 -2.8320 3.7760 2.8320 -2.8320 -2.1240 2.8320 2.1240k4=Bar2D2Node_Stiffness (E,A,x4,y4,x3,y3,alpha1) k4 = 73750000 0 -73750000 0 0 0 0 0 -73750000 0 73750000 0 0 0 0 0用到的函数直接复制到下面的目录：C:Documents and SettingslidejunMy DocumentsMATLABfunction k=B

9、ar2D2Node_Stiffness(E,A,x1,y1,x2,y2,alpha)%该函数计算单元的刚度矩阵%输入弹性模量E，横截面积A%输入第一个节点坐标（x1,y1），第二个节点坐标（x2,y2），角度alpha（单位是度）%输出单元刚度矩阵k(4X4)。%-L=sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);x=alpha*pi/180;C=cos(x);S=sin(x);k=E*A/L*C*C C*S -C*C -C*S; C*S S*S -C*S -S*S; -C*C -C*S C*C C*S; -C*S -S*S C*S S*S;（3） 建立整体刚度方程由

10、于该结构共有4个节点，因此，设置结构总的刚度矩阵为KK(8×8)，先对KK清零，然后四次调用函数Bar2D2Node _Assembly进行刚度矩阵的组装。相关的计算流程如下。KK=zeros(8,8);KK=Bar2D2Node_Assembly (KK,k1,1,2);KK=Bar2D2Node_Assembly (KK,k2,2,3);KK=Bar2D2Node_Assembly (KK,k3,1,3);KK=Bar2D2Node_Assembly (KK,k4,4,3)KK= 1.0e+008 * 1.1151 0.2832 -0.7375 0 -0.3776 -0.2832

11、 0 0 0.2832 0.2124 0 0 -0.2832 -0.2124 0 0 -0.7375 0 0.7375 0.0000 -0.0000 -0.0000 0 0 0 0 0.0000 0.9833 -0.0000 -0.9833 0 0 -0.3776 -0.2832 -0.0000 -0.0000 1.1151 0.2832 -0.7375 0 -0.2832 -0.2124 -0.0000 -0.9833 0.2832 1.1957 0 0 0 0 0 0 -0.7375 0 0.7375 0 0 0 0 0 0 0 0 0function z = Bar2D2Node_Ass

12、embly(KK,k,i,j)%该函数进行单元刚度矩阵的组装%输入单元刚度矩阵k，单元的节点编号i、j%输出整体刚度矩阵KK%-DOF(1)=2*i-1;DOF(2)=2*i;DOF(3)=2*j-1;DOF(4)=2*j;for n1=1:4 for n2=1:4 KK(DOF(n1),DOF(n2)= KK(DOF(n1),DOF(n2)+k(n1,n2); endendz=KK;（4） 边界条件的处理及刚度方程求解 由图3-8可以看出，节点1的位移将为零，即, ，节点2的位移，节点4的，。节点载荷10N。采用高斯消去法进行求解，注意：MATLAB中的反斜线符号“”就是采用高斯消去法。该结

13、构的节点位移为：而节点力为：其中，为节点1处沿x和y方向的支反力，为节点2处y方向的支反力，为节点4处沿x和y方向的支反力。相关的计算流程如下。k=KK(3,5,6,3,5,6) k =1.0e+008 * 0.7375 -0.0000 -0.0000 -0.0000 1.1151 0.2832 -0.0000 0.2832 1.1957p=20000;0;-25000;u=kpu = 1.0e-003 * 0.2712 0.0565 -0.2225 这里将列排成了一行，以节省篇幅排量 111111111111111111111111111111111111111111111111111111

14、111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111 由此可以看出，所求得的结果，则所有节点位移为 (3-75)与前面通过数学推导得到的结果相同，见式(3-72)。（5）支反力的计算在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。将整体的位移列阵

15、q（采用国际单位）代回原整体刚度方程，计算出所有的节点力P，按上面的对应关系就可以找到对应的支反力。相关的计算流程如下。q=0 0 0.0002712 0 0.0000565 -0.0002225 0 0'q = 1.0e-003 * 0 0 0.2712 0 0.0565 -0.2225 0 0 这里将列排成了一行，以节省篇幅排量 11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111

16、1111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111P=KK*qP = 1.0e+004 * -1.5833 0.3126 2.0001 2.1879 -0.0001 -2.5005 -0.4167 0 这里将列排成了行 按对应关系，可以得到对应的支反力为 (3-76)（6）各单元的应力计算先从整体位移列阵q中提取出单元的位移列阵，然后，调用计算单元应力的函数Bar2D2Node_ElementStress，就可

17、以得到各个单元的应力分量。当然也可以调用上面的Bar2D2Node_ElementForces(E,A,x1,y1,x2,y2,alpha,u)函数来计算单元的集中力，然后除以面积求得单元应力。相关的计算流程如下。>>u1=q(1);q(2);q(3);q(4)u1 = 1.0e-003 * 0 0 0.2712 0 这里将列排成了一行，以节省篇幅排量 12121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212

18、1212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212>> stress1=Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha1,u1) stress1 =2.0001e+008>>u2=q(3);q(4);q(5);q(6)u2 = 1.0e-003 * 0.2712 0 0.0565 -0.2225 这里将列排成了一行，以节省篇幅排量 12

19、1212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212>> stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,alpha2,

20、u2) stress2 = -2.1879e+008>>u3=q(1);q(2);q(5);q(6)u3 = 1.0e-003 * 0 0 0.0565 -0.2225 这里将列排成了一行，以节省篇幅排量 12121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121

21、2121212121212121212121212121212121212121212121212121212>> stress3=Bar2D2Node_Stress(E,x1,y1,x3,y3,alpha3,u3) stress3 = -52097000>>u4=q(7);q(8);q(5);q(6)u4 = 1.0e-003 * 0 0 0.0565 -0.2225 这里将列排成了一行，以节省篇幅排量 12121212121212121212121212121212121212121212121212121212121212121212121212121212121

22、2121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212121212>> stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,alpha1,u4) stress4 = 41668750function stress= Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)%该

23、函数计算单元的应力%输入弹性模量E，第一个节点坐标（x1,y1），第二个节点坐标（x2,y2）%输入角度alpha（单位是度）和单位节点位移矢量u%返回单元应力标量stress %- L=sqrt(x2-x1)*(x2-x1)+(y2-y1)*(y2-y1);x=alpha*pi/180;C=cos(x);S=sin(x);stress=E/L*-C -S C S*u;可以看出：计算得到的单元1的应力为；单元2的应力为，单元3的应力为，单元4的应力为。与前面通过数学推导得到的结果相同。1.5. MATLAB直接使用FEM求解Excel中的简单温度场分布>>pdetool具体操作步骤

24、：1. 画出正方形，单击后设定其大小及位置2. 进入Boundary模式，设定边界条件boundary condition1）边界条件是先显示边的编号，选中该边界在双击即可对每个边界输入参数。2) 固定边界条件选择R表示真实的值，h表示系数，u就是温度分别设定最上边的边的温度r=100, 左边r=75, 右边r=50, 底边r=03. 进入PDE菜单，指定PDE参数（PDE Specification）静态的温度场的求解是椭圆型Possion方程而我们对比下来，可以知道是f=0, a=0, c=1即：4. 网格划分(Mesh, Initialize Mesh)5. 求解Solve( Solve

25、 PDE)1.6. 有限元的线性方程组的求解一）直接消元法1. 最简单的当然是高斯消元法，为了保证精度，可以取主元2. LU三角分解，对应的也就是高斯消元法二）迭代法1.雅可比迭代法： （i=1,2,3n）先给出一组的初始值（i=1,2,3n）然而按照公式进行迭代就可以得到，采用Excel求解得到的迭代是：Yx0x15-310200.4-311-311401.2727271-311-31601.45454501-361702.8333335-31020.40.872727-311-31141.2727271.5209371-311-3161.4545452.53801701-36172.833

26、3333.3484855-31020.8727270.804959-311-31141.5209371.8985221-311-3162.5380172.70323101-36173.3484853.8488522.高斯-赛德尔迭代雅可比迭代中在计算的过程中，总是使用上一次的计算值，赛德尔迭代中，每次算的时候尽量使用最新值 （i=1,2,3n）Yx0x15-310200.4-311-311401.3818181-311-31601.79504101-361703.5005515-31020.40.870083-311-31141.3818181.6813471-311-3161.7950412

27、.78869201-36173.5005513.9474555-31020.8700830.85107-311-31141.6813471.906531-311-3162.7886922.97371701-36173.9474554.0024375-31020.851070.949175-311-31141.906531.9787491-311-3162.9737172.99948901-36174.0024374.003286比较上面的值，可以看出高斯-赛德尔迭代收敛速度比雅克比迭代快太多了。2. 有限元建模2.1. ANSYS10.0的起始菜单的说明：topological Opt:拓扑优

28、化Design Opt: 设计优化Radiation Opt: 辐射优化Run-Time Stats: 运行状态ansys软件的probdesign (概率设计)Chapter 6. Reduced Order Modeling (ROM):This chapter describes a solution method for efficiently solving coupled-field problems involving flexible structures. This reduced order modeling (ROM) method is based on a modal

29、 representation of the structural response. The deformed structural domain is described by a factored sum of the mode shapes (eigenvectors). The resulting ROM is essentially an analytical expression for the response of a system to any arbitrary excitation.This methodology has been implemented for co

30、upled electrostatic-structural analysis and is applicable to micro-electromechanical systems (MEMS).The solver tool enables essential speed up for two reasons:A few eigenmodes accurately represents the dynamic behavior of most structures. This is particularly true for micro-electromechanical systems

31、 (MEMS).Modal representations of electrostatic-structural domains are very efficient because just one equation per mode and one equation per conductor are necessary to describe the coupled domain system entirely.This modal method can be applied to nonlinear systems. Geometrical nonlinearities, such

32、as stress stiffening, can be taken into account if the modal stiffness is computed from the second derivatives of the strain energy with respect to modal coordinates. Capacitance stroke functions provide nonlinear coupling between eigenmodes and the electrical quantities if stroke is understood to b

33、e modal amplitude.For more information on this method, see Reduced Order Modeling of Coupled Domains in the ANSYS, Inc. Theory Reference. The process involves the following distinct steps.2.2. ANSYS仿真中常用的单元：# 最常见的solid45单元SOLID45 Element DescriptionSOLID45 is used for the 3-D modeling of solid struc

34、tures. The element is defined by eight nodes having three degrees of freedom at each node: translations in the nodal x, y, and z directions.The element has plasticity, creep, swelling, stress stiffening, large deflection, and large strain capabilities. A reduced integration option with hourglass con

35、trol is available. See SOLID45 in the ANSYS, Inc. Theory Reference for more details about this element. A similar element with anisotropic properties is SOLID64. A higher-order version of the SOLID45 element is SOLID95.SOLID45 Input DataThe geometry, node locations, and the coordinate system for thi

36、s element are shown in Figure 45.1: "SOLID45 Geometry". The element is defined by eight nodes and the orthotropic(各向异性) material properties. Orthotropic material directions correspond to the element coordinate directions. The element coordinate system orientation is as described in Coordin

37、ate Systems.Element loads are described in Node and Element Loads. Pressures may be input as surface loads on the element faces as shown by the circled numbers on Figure 45.1: "SOLID45 Geometry". Positive pressures act into the element. Temperatures and fluences may be input as element bod

38、y loads at the nodes. The node I temperature T(I) defaults to TUNIF. If all other temperatures are unspecified, they default to T(I). For any other input temperature pattern, unspecified temperatures default to TUNIF. Similar defaults occurs for fluence except that zero is used instead of TUNIF.KEYO

39、PT(1) is used to include or suppress the extra displacement shapes. KEYOPT(5) and KEYOPT(6) provide various element printout options (see Element Solution).This element also supports uniform reduced (1 point) integration with hourglass control when KEYOPT(2) = 1. Using uniform reduced integration pr

40、ovides the following advantages when running a nonlinear analysis: Less CPU time is required for element stiffness formation and stress/strain calculations to achieve a comparable accuracy to the FULL integration option.The length of the element history saved record (.ESAV and .OSAV) is about 1/7th

41、as much as when the full integration (2 X 2 X 2) is used for the same number of elements.Nonlinear convergence characteristic of the option is generally far superior to the default full integration with extra displacement shape; that is, KEYOPT(1) = 0, KEYOPT(2) = 0.The analysis will not suffer from

42、 volumetric locking which can be caused by plasticity or other incompressible material properties.An analysis using uniform reduced integration can have the following disadvantages: The analysis is not as accurate as the full integration method, which is apparent in the linear analysis for the same

43、mesh.The analysis cannot capture the bending behavior with a single layer of elements; for example, in the case of a fixed-end cantilever with a lateral point load, modeled by one layer of elements laterally. Instead, four elements are usually recommended.When the uniform reduced integration option

44、is used (KEYOPT(2) = 1 - this option is the same as SOLID185 with KEYOPT(2) = 1), you can check the accuracy of the solution by comparing the total energy (SENE label in ETABLE) and the artificial energy (AENE label in ETABLE) introduced by hourglass control. If the ratio of artificial energy to tot

45、al energy is less than 5%, the solution is generally acceptable. If the ratio exceeds 5%, refine the mesh. The total energy and artificial energy can also be monitored by using the OUTPR,VENG command in the solution phase. For more details, see the ANSYS, Inc. Theory Reference.You can apply an initi

46、al stress state to this element through the ISTRESS or ISFILE command. For more information, see Initial Stress Loading in the ANSYS Basic Analysis Guide. Alternately, you can set KEYOPT(9) = 1 to read initial stresses from the user subroutine USTRESS. For details on user subroutines, see the Guide

47、to ANSYS User Programmable Features.# 最常见的PlANE82单元PLANE82PLANE82 Element DescriptionPLANE82 is a higher order version of the 2-D, four-node element (PLANE42). It provides more accurate results for mixed (quadrilateral-triangular) automatic meshes and can tolerate irregular shapes without as much lo

48、ss of accuracy. The 8-node elements have compatible displacement shapes and are well suited to model curved boundaries.The 8-node element is defined by eight nodes having two degrees of freedom at each node: translations in the nodal x and y directions. The element may be used as a plane element or

49、as an axisymmetric element. The element has plasticity, creep, swelling, stress stiffening, large deflection, and large strain capabilities. Various printout options are also available. See PLANE82 in the ANSYS, Inc. Theory Reference for more details about this element. See PLANE83 for a description

50、 of an axisymmetric element which accepts nonaxisymmetric loading.The geometry, node locations, and the coordinate system for this element are shown in Figure 82.1: "PLANE82 Geometry".A triangular-shaped element may be formed by defining the same node number for nodes K, L and O. A similar

51、, but 6-node, triangular element is PLANE2. Besides the nodes, the element input data includes a thickness (TK) (for the plane stress option only) and the orthotropic material properties. Orthotropic material directions correspond to the element coordinate directions. The element coordinate system orientation is as described in Coordinate Systems.Element loads are described in Node and Element Loads. Pressures may be input as surface loads on the element faces as shown by the circled numbers on Figure 82.1: "PLANE82 Geometry". Positive pressures act into the element. Temperatures an