ANSYS的仿真及教学日记
1. 有限元的基础知识 (2)
1.1. 函数逼近的思想 (2)
1.2. 差分法求解偏微分方程 (4)
1.3. 差分法求解偏微分方程............................................................. 错误！未定义书签。
1.4. 等参变换、变分法、单元刚度矩阵相关 (4)
1.5. MATLAB有限元编程求解力学问题 (7)
1.6. MATLAB直接使用FEM求解Excel中的简单温度场分布 (12)
1.7. 有限元的线性方程组的求解 (13)
2. 有限元建模 (16)
2.1. ANSYS10.0的起始菜单的说明： (16)
2.2. ANSYS仿真中常用的单元： (17)
3. 静力学分析 (24)
3.1. 壳体结构的有限元分析 (24)
3.2. 几个教学案例的分析要点： (24)
3.3. 结构分析线性方程求解的基本算法 (28)
1.有限元的基础知识
1.1.函数逼近的思想
(a) 基于全域[x0,x L]的函数逼近（全域高阶插值或拟合的思想）
(b) 基于子域[x i,x i+1]的函数逼近（子域低阶分段插值或分段拟合的思想-有限元的思想）
1.2. 差分法求解偏微分方程
Excel 演示：
上边（100℃,左为75℃，右为50℃，底为0℃）的温度场
1.3. 等参变换、变分法、单元刚度矩阵相关
#.直接刚度矩阵建立与求解 1) 什么是应力、应变、位移
2）什么是插值函数？什么是高次单元？什么是低次单元？ 3）直接刚度法？
4）有限元法求解的是什么方程？
#.刚度与刚度矩阵
弹簧的变形与力的关系
F l k =??
k 是弹簧系数，是弹簧的固有特性，和截面积、长度、形状有关系
变成矩阵的形式： y x A =? 也就是
e e e F q K =?
对于杆件：
()T
e u u q 21=，节点位移
()T
e F F F 21=，节点外力
该节点的内力为：
????????????????=??????2121---u u l EA l
EA l
EA l EA F F I I （正负号是考虑到方向） 而内力和外力的平衡有11-I F F =，22I F F =, 因此该方程
??????=????????????????2121--F F u u l EA l
EA l EA l EA
???????
???=l EA l
EA l EA l EA K e --就是刚度矩阵
#.等参变换
由于子单元的形状复杂多变，采用等参变换后的坐标的刚度矩阵除了一些非常简单的情况，基本上不能进行显式积分，而需要进行数值积分，如采用Newton-Cotes 积分等等，基于对积分方程的各种假设，比如使用二次的线性来代替，或者三次线性函数代替，或者其他假设的函数等等，这样的单元刚度矩阵有各种各样的方法了。
3. 单元的各个参量的插值函数
这是一个自由度的情况，该单元位移场为:
() +++=2210x x x U ααα
取首两项作低阶插值：
()x x U 10αα+=
节点条件为：
()10u x u x == ()2u x u l x ==
带入得到
10u =α，()
l
u u 121-=
α
带入得到：
()
()e x q x N x l
x u l x x l
u u u u =+??? ??=+
=1121-1- N 叫做形状函数矩阵，具体为
()??
? ??=l x l
x
x N -1 q 是节点位移矩阵
()T
e u u q 21=
#. 整体坐标系和局部坐标系
#. 能量原理（虚位移原理、最小势能原理）、变分法和有限元的建立 物理方程就是微分方程加上边界条件，对微分方程的原方程的求解是非常困难，甚至不可能的，变分法就是先假设一个泛函（这个泛函是容易处理的标准函数，如线性函数等），泛函的求导就是这个微分方程，通过泛函的求极值来满足边界条件，最后还是得到单元刚度矩阵。
通过胡克定律和能量方程等可以求出中性层的挠度方程为
())2-(240?,3340
x l lx x EI
p y
x v +== （p 为压强）
最小势能原理得到（假设一个正弦函数为基函数）： 应变能为： W=
1.4. MATLAB 有限元编程求解力学问题
四杆桁架结构的有限元分析(Bar2D2Node)
如图所示的结构，各个杆的弹性模量和横截面积都为4229.510/E N mm =?, 2100A mm =。试基于MATLAB 平台求解该结构的节点位移、单元应力以及支反力。
图四杆桁架结构
解答：对该问题进行有限元分析的过程如下。
（1）结构的离散化与编号
对该结构进行自然离散，节点编号和单元编号如图。
（2）计算各单元的刚度矩阵(基于国际标准单位)
建立一个工作目录，将所编制的用于平面桁架单元分析的4个MATLAB函数放置于该工作目录中，分别以各自函数的名称给出文件名，即：Bar2D2Node_Stiffness，Bar2D2Node_Assembly，Bar2D2Node_Stress，Bar2D2Node_Forces。然后启动MATLAB，将工作目录设置到已建立的目录中，在MATLAB环境中，输入弹性模量E、横截面积A，各点坐标x1,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 0
0 0 0 0
-73750000 0 73750000 0
0 0 0 0
k2=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.8333
k3=Bar2D2Node_Stiffness (E,A,x1,y1,x3,y3,alpha3)
k3 = 1.0e+007 *
3.7760 2.8320 -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.1240
k4=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 Settings\lidejun\My Documents\MATLAB
function k=Bar2D2Node_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）建立整体刚度方程
由于该结构共有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 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 0
function z = Bar2D2Node_Assembly(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); end end z=KK;
（4） 边界条件的处理及刚度方程求解
由图3-8可以看出，节点1的位移将为零，即10u =, 10v =，节点2的位移20v =，节点4的40u =，40v =。节点载荷3F ＝10N 。采用高斯消去法进行求解，注意：MATLAB 中的反斜线符号“\”就是采用高斯消去法。
该结构的节点位移为：[
]T
v u v u v u v u 44
3322
11=q
而节点力为：[
]
T
y x y y x R R R R R 444
2
4
1
1
105.20102?-?=+=F R P
其中，11(,)x y R R 为节点1处沿x 和y 方向的支反力，2y R 为节点2处y 方向的支反力，
44(,)x y R R 为节点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.1957
p=[20000;0;-25000]; u=k\p
u = 1.0e-003 *
0.2712 0.0565 -0.2225 [这里将列排成了一行，以节省篇幅]
由此可以看出，所求得的结果2330.271 2,0.056 5,0.222 5u mm u mm v mm ===-，则所有节点位移为
[]000.271 200.056 50.222 500T
mm =-q (3-75)
与前面通过数学推导得到的结果相同，见式(3-72)。
（5）支反力的计算
在得到整个结构的节点位移后，由原整体刚度方程就可以计算出对应的支反力。将整体的位移列阵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 [这里将列排成了一行，以节省篇幅] P=KK*q
P = 1.0e+004 *
-1.5833 0.3126 2.0001 2.1879 -0.0001 -2.5005 -0.4167 0 [这里将列排成了行]
按对应关系，可以得到对应的支反力为
(3-76)
（6）各单元的应力计算
先从整体位移列阵q 中提取出单元的位移列阵，然后，调用计算单元应力的函数Bar2D2Node_ElementStress ，就可以得到各个单元的应力分量。当然也可以调用上面的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 [这里将列排成了一行，以节省篇幅] >> 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 [这里将列排成了一行，以节省篇幅]
>> stress2=Bar2D2Node_Stress(E,x2,y2,x3,y3,alpha2,u2) stress2 = -2.1879e+008
>>u3=[q(1);q(2);q(5);q(6)] u3 = 1.0e-003 *
0 0 0.0565 -0.2225 [这里将列排成了一行，以节省篇幅] >> 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 [这里将列排成了一行，以节省篇幅] >> stress4=Bar2D2Node_Stress(E,x4,y4,x3,y3,alpha1,u4) stress4 = 41668750
function stress= Bar2D2Node_Stress(E,x1,y1,x2,y2,alpha,u)
%该函数计算单元的应力
%输入弹性模量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的应力为82.000 110Pa e σ=?； 单元2的应力为82.187 910Pa σ=-?， 单元3的应力为75.209 710Pa σ=-?，
单元4的应力为74.16710Pa σ=?。与前面通过数学推导得到的结果相同。
1.5. MATLAB 直接使用FEM 求解Excel 中的简单温度场分布
>>pdetool
具体操作步骤：
1. 画出正方形，单击后设定其大小及位置
2. 进入Boundary 模式，设定边界条件boundary condition
1）边界条件是先显示边的编号，选中该边界在双击即可对每个边界输入参数。 2) 固定边界条件选择
R 表示真实的值，h 表示系数，u 就是温度
分别设定最上边的边的温度r=100, 左边r=75, 右边r=50, 底边r=0
3. 进入PDE 菜单，指定PDE 参数（PDE Specification ）
静态的温度场的求解是椭圆型Possion 方程
()()0-=???
?
??????+??? ??????=
???=+???y u c y u x u c x u u c f au u c 而我们
0222
2=??+??y
T
x T
对比下来，可以知道是f=0, a=0, c=1 即：
4. 网格划分(Mesh, Initialize Mesh)
5. 求解Solve( Solve PDE)
1.6. 有限元的线性方程组的求解
一）直接消元法
1. 最简单的当然是高斯消元法，为了保证精度，可以取主元
2. LU 三角分解，对应的也就是高斯消元法
二）迭代法
1.雅可比迭代法：
????
? ??=
∑≠=n i j j j ij i ii
i x a b a x 1-1 （i=1,2,3…n ）
先给出一组i x 的初始值0i x （i=1,2,3…n ）然而按照公式进行迭代就可以得到
b Ax =
?
?
??
?
?
???
???=63-103-113-113-113-013-5A ，????????????=1716142b
2.高斯-赛德尔迭代
雅可比迭代中在计算1+k i x 的过程中，总是使用上一次的计算值k i x ，赛德尔迭代中，每次算的时候尽量使用最新值
???
?
??=∑∑+==++n i j k j ij i j k j ij i ii
k i
x a x a b a x
11-11
1--1 （i=1,2,3…n ）
比较上面的值，可以看出高斯-赛德尔迭代收敛速度比雅克比迭代快太多了。
2.有限元建模
2.1.ANSYS10.0的起始菜单的说明：
topological Opt:拓扑优化
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 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 coupled 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 (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 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 be 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 Description
SOLID45 is used for the 3-D modeling of solid structures. 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 control 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 Data
The geometry, node locations, and the coordinate system for this 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 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 45.1: "SOLID45 Geometry". Positive pressures act into the element. Temperatures and fluences may be input as element body 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.
KEYOPT(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 provides 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 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 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 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 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 total 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 initial 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 to ANSYS User Programmable Features.
# 最常见的PlANE82单元
PLANE82
PLANE82 Element Description
PLANE82 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 loss 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 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 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, 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