基于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)

k2 =

+006 *

Columns 1 through 5

0 0

Column 6

3.集成整体刚度矩阵 8*8零矩阵

K =

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 >> K=LinearTriangleAssemble(K,k1,1,3,4)

K =

+006 *

Columns 1 through 5

0 0 0 0

0 0 0

0 0 0 0 0

0 0 0 0 0 0 0 0

Columns 6 through 8

0 0 0 0 0 0 0

>> K=LinearTriangleAssemble(K,k1,1,2,3) K =

+007 *

0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0

0 0 0 0

4.引入边界条件.用上一步得到的整体刚度矩阵，可以得到该结构的方程组如下形式

??????????

??????????????=??????????????????

??????????

????????????????????????4y 4X 3y 3X 2y 2X 1y 1X 4y 4X 3y 3X 2y 2X 1y 1X 6F F F F F F F F U U U U U U U U 6.2740 1.8750- 0.5048- 0.8654 0 0 5.7692- 1.0096 1.8750- 3.4615 1.0096 1.4423- 0 0 0.8654 2.0192- 0.5048- 1.0096 6.2740 0 5.7692- 0.8654 0

1.8750- 0.8654 1.4423- 0 3.4615 1.0096

2.0192- 1.8750- 0 0 0 5.7692- 1.0096 6.2740 1.8750- 0.5048- 0.8654 0 0 0.8654 2.0192- 1.8750-

3.4615 1.0096 1.4423- 5.7692- 0.8654 0 1.8750- 0.5048- 1.0096 6.2740 0 1.0096 2.0192- 1.8750- 0 0.8654 1.4423- 0

3.461510 本题的边界条件：

04411====y x y x U U U U

0,375.9,0,375.93322====y x y x F F F F

将边界条件带入，得到：

??????

???

?????

???????

?????=???????????????

???

?????

???????????????????????

??????4y 4X 1y 1X 3y 3X 2y 2X 6F F 09.37509.375F F 0 0 U U U U 0 0 6.2740 1.8750- 0.5048- 0.8654 0 0 5.7692- 1.0096 1.8750- 3.4615 1.0096 1.4423- 0 0 0.8654 2.0192- 0.5048- 1.0096 6.2740 0 5.7692- 0.8654 0

1.8750- 0.8654 1.4423- 0 3.4615 1.0096

2.0192- 1.8750- 0 0 0 5.7692- 1.0096 6.2740 1.8750- 0.5048- 0.8654 0 0 0.8654 2.0192- 1.8750-

3.4615 1.0096 1.4423- 5.7692- 0.8654 0 1.8750- 0.5048- 1.0096 6.2740 0 1.0096 2.0192- 1.8750- 0 0.8654 1.4423- 0

3.461510

5.解方程

分解上述方程组，提取总体刚度矩阵K 的第3-6行的第3-6列作为子矩阵

????

???=?????????????????????????? 09.37509.375 U U U U 6.2740 0 5.7692- 0.8654 0 3.4615 1.0096 2.0192- 5.7692- 1.0096 6.2740 1.8750- 0.8654 2.0192- 1.8750-

3.4615103y 3X 2y 2X 6 Matlab 命令

>> k=K(3:6,3:6) k =

+006 *

0 0

>> f=[;0;;0] f = 0 0 >> u=k\f u =

现在可以清楚的看出，节点2的水平位移和垂直位移分别是0.7111m 和0.1115m 。节点3的水平位移和垂直位移分别是0.6531m 和0.0045m 。 6.后处理

用matlab 命令求出节点1和节点4的支反力以及每个单元的应力。首先建立总体节点位移矢量U ，

U=[0;0;u;0;0] U = *

0 0 0 0 >> F=K*U F =

由以上知，节点1的水平反力和垂直反力分别是（指向左边）和（作用力方向向下），节点4的水平反力和垂直反力分别是（指向左边）和（作用力方向向下）.满足力平衡条件。接着，建立单元节点位移矢量21u 和u ，然后调用matlab 命令LinearTriangleElementStresses 计算单元应力sigma1和sigma2

>> u1=[U(1);U(2);U(5);U(6);U(7);U(8)]

u1 =

>> u2=[U(1);U(2);U(3);U(4);U(5);U(6)]

u2 =

>> sigma1=LinearTriangleElementStresses(E,NU,,0,0,,,0,,1,u1) sigma1 =

+003 *

>> sigma2=LinearTriangleElementStresses(E,NU,,0,0,,0,,,1,u2) sigma2 =

+003 *

由以上可知，单元1的应力)(0144.3拉应力MPa x =σ，)(9043.0y 拉应力MPa =σ，

)(0072.0y 正值MPa x =τ 。单元2的应力是

)(9856.2拉应力MPa x =σ)(0036.0y 压应力MPa =σ)(0072.0y 负值MPa x =τ。显然，在x

方向的应力（拉应力）接近于正确的值3MPa （拉应力）。接着调用

LinearTriangleElementStresses 函数计算每个单元的主应力和主应力方向角。 >> s1= LinearTriangleElementPStresses(sigma1) s1 = +003 *

>> s2= LinearTriangleElementPStresses(sigma2) s2 = +003 *

)(0144.31拉应力MPa =σ，)(9043.02拉应力MPa =σ主应力方向角ο2.0=p θ )(M 9856.21拉应力Pa =σ，)(0036.02压应力MPa =σ，ο1.0-=p θ

例 2.考虑如图所示的由均匀分布载荷和集中载荷作用的薄平板结构。将平板离散化成12个线性三角单元，如图4所示。假定E=210GPa

，

t=0.025m

，

w=100kN/m

和

P=。

1.离散化

2.写出单元刚度矩阵

>> E=201e6;

>> NU=;

>> t=;

>> k1= LinearTriangleElementStiffness(E,NU,t,0,,,,,,1); >> k2= LinearTriangleElementStiffness(E,NU,t,0,,0,,,,1); >> k3= LinearTriangleElementStiffness(E,NU,t,,,,,,,1); >> k4= LinearTriangleElementStiffness(E,NU,t,,,0,,,,1); >> k5= LinearTriangleElementStiffness(E,NU,t,0,,,,,,1); >> k6= LinearTriangleElementStiffness(E,NU,t,0,,0,0,,,1); >> k7= LinearTriangleElementStiffness(E,NU,t,,,,,,0,1); >> k8= LinearTriangleElementStiffness(E,NU,t,,,0,0,,0,1); >> k9= LinearTriangleElementStiffness(E,NU,t,025,,,0,,,1); >> k10= LinearTriangleElementStiffness(E,NU,t,,,,,,,1); >> k11= LinearTriangleElementStiffness(E,NU,t,,0,,0,,,1); >> k12= LinearTriangleElementStiffness(E,NU,t,,,,0,,,1)

k1 =

+006 *

3.集成整体刚度矩阵：

>>K=zero(22,22);

>>K=LinearTriangleAssemble(K,k1,1,3,2);

>>K=LinearTriangleAssemble(K,k2,1,4,3);

>>K=LinearTriangleAssemble(K,k3,3,5,2);

>>K=LinearTriangleAssemble(K,k4,3,4,5);

>>K=LinearTriangleAssemble(K,k5,4,6,5);

>>K=LinearTriangleAssemble(K,k6,4,7,6);

>>K=LinearTriangleAssemble(K,k7,5,6,8);

>>K=LinearTriangleAssemble(K,k8,6,7,8);

>>K=LinearTriangleAssemble(K,k9,5,8,9);

>>K=LinearTriangleAssemble(K,k10,5,9,10);

>>K=LinearTriangleAssemble(K,k11,8,11,9);

>>K=LinearTriangleAssemble(K,k12,9,11,10)

运行得

+008 *

Columns 1 through 7

0 0

0 0 0

0 0

0 0 0

0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Columns 8 through 14

0 0 0 0 0 0

0 0 0 0 0

0 0 0 0

0 0

0 0 0

0 0

0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0

Columns 15 through 21

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0 0 0

0 0

Column 22

4.引入边界条件：

U1x= U1y= U4x= U4y=U7x=U7y=0

F2x= F2y= F3x= F3y=F6x=F6y=F8x= F8y= F9x= F9y=F10x=F10y= F11x= F11y= 0

F5x= 0,F5y=

5.解方程：

>>k=[K(3:6,3:6),K(3:6,9:12),K(3:6,15:22);K(9:12,3:6),K(9:12,9:12),K(9:12,15:22) ;K(15:22,3:6),K(15:22,9:12) ,K(15:22,15:22)];

K =

+008 *

Columns 1 through 8

0 0

0 0 0

0 0 0 0 0 0

0 0 0 0 0 0 0 0

Columns 9 through 16

0 0 0 0 0 0 0 0

0 0 0 0 0 0

0 0 0

0 0

0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Columns 17 through 22

0 0 0 0 0 0

0 0

0 0 0 0 0 0

0 0

>>f=[0;0;0;0;0;;0;0;0;0;0;0;0;0;0;0];

>>u=k\f

*[ ]T

6.后处理

>>U=[0;0;u(1:4);0;0;u(5:8);0;0;u(9:16)];

>>F=K*U

[ 0 0 ]T

三.小结

通过这次练习，对有限元结合MATLAB软件解决平面应力/应变问题的方法和过程逐渐清晰，特别是在应用MATLAB软件的过程中学到了很多东西。同时，通过此次课程设计，我对以前学过的有限元理论课程知识也有了进一步了解。一开始还没有头绪，耗了半天也没正处东西来，后来吧课本拿来仔细看模仿学习，逐渐对软件有所熟悉，对解题思路更加清晰，我深刻的认识到自己动手自己找方法的重要性。以后要加强这方面的能力。

matlab有限元分析实例

MATLAB: MATLAB是美国MathWorks公司出品的商业数学软件，用于数据分析、无线通信、深度学习、图像处理与计算机视觉、信号处理、量化金融与风险管理、机器人，控制系统等领域。 MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室），软件主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等。MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。 MATLAB有限元分析与应用:

《MATLAB有限元分析与应用》是2004年4月清华大学出版社出版的图书，作者是卡坦，译者是韩来彬。内容简介：《MATLAB有限元分析与应用》特别强调对MATLAB的交互应用，书中的每个示例都以交互的方式求解，使读者很容易就能把MATLAB用于有限分析和应用。另外，《MATLAB有限元分析与应用》还提供了大量免费资源。《MATLAB有限元分析与应用》采用当今在工程和工程教育方面非常流行的数学软件MATLAB来进行有限元的分析和应用。《MATLAB有限元分析与应用》由简单到复杂，循序渐进地介绍了各种有限元及其分析与应用方法。书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题，具有很高的工程应用价值。

局部应力应变分析法

1.局部应力应变分析法、名义应力疲劳设计法、疲劳可靠性设计法、损伤容限设计法 2.磨损、腐蚀、断裂 3.交变应力水平低、脆性断裂、损伤积累过程、断口在宏观和微观上有特征 4.表面应力水平比内部高、表面晶体束缚少，易发生滑移、表面易发生环境介质腐蚀、表面的加工痕迹或划痕会降低零件疲劳强度 5.材料在循环应力、应变作用下，某点或某些点发生局部永久性结构变形，在经过一定循环次数后产生裂纹或发生断裂的过程。 6.外加应力水平和标准试样疲劳寿命之间关系的曲线 7.疲劳寿命无穷大时的中值疲劳强度 8.在各级应力水平下的疲劳寿命分布曲线上可靠度相等的点连成曲线就能得到给定可靠度的一组SN曲线 9.理论应力:局部应力与名义应力的比值Kt=6t/6n 10.在应力集中和终加工相同的情况下，尺寸为d的零件的极限寿命与标准直径试样的极限寿命的比值 11.史密斯图、海夫图、等寿命图（相同寿命时在不同应力下的疲劳极限间关系的线图） 12.线性积累损伤理论： 13.载荷随时间变化的历程应力随时间变化的历程 14.零件的疲劳破损都是从应变集中部位最大局部应变处开始的裂纹萌生以前，一般都会产生塑性变形塑性变形是裂纹萌生和扩展的先决条件零件的疲劳强度和寿命由应变集中部位的最大局部应力应变决定 15参数应力（名义应力）应变（局部应变）特征应力疲劳应变疲劳范围104-105-5*106 103-104-105 寿命总寿命裂纹形成寿命曲线SN曲线古德曼曲线EN曲线，循环应力应变曲线变形弹性变形应力应变成正比塑性变形较大 16真实应力 17材料在循环载荷作用下的应力应变响应循环应力应变曲线 18循环硬化：应力幅6a为常数，应变幅Ea随着循环次数增加而减少，最后趋于稳定循环软化：应变幅Ea为常数，应力幅6a随着循环次数增加而逐渐减少 19.漫森四点：应变寿命曲线的弹性线上取2点，塑性线上取2点，通用斜率法 20.雨流法：Y方向为时间，X方向为应力大小 21.在循环加载作用下应力应变响应称为循环应力应变曲线在循环加载作用下应力应变轨迹线称为应力应变迟滞回线件加载拉伸到A卸载到O加载压缩到B加载拉伸到C（与A重合）形成的环线 22.损伤容限设计：以断裂力学理论为基础以无损检测技术和断裂韧性与疲劳裂纹扩展速率的测定技术为手段以有初始缺陷的寿命估算为中心以断裂控制为保障确保零件在使用期内能够安全使用的一种疲劳计算方法 23.应力强度因子：K是度量裂纹端部应力场强弱程度的一个参数 24.断裂韧度：应力强度因子的临界值，发生脆断时的应力强度因子。 25.性能、可靠性（规定条件规定时间完成规定功能）、维修性指标（规定条件时间程序方法恢复到规定状态） 26.广义可靠性=狭义可靠性（不可维修产品的可靠性）+可维修性 27.故障和失效（产品不能完成其规定功能的状态） 28.可靠度（规定条件时间完成规定功能的概率）

Matlab有限元分析操作基础

Matlab 有限元分析20140226 为了用Matlab 进行有限元分析，首先要学会Matlab 基本操作，还要学会使用Matlab 进行有限元分析的基本操作。 1. 复习：上节课分析了弹簧系统 x 推导了系统刚度矩阵 11221 21200k k k k k k k k -????-????--+??

2. Matlab有限元分析的基本操作（1）单元划分（选择何种单元，分成多少个单元，标号）（2）构造单元刚度矩阵（列出…）（3）组装系统刚度矩阵（集成整体刚度矩阵）（4）引入边界条件（消除冗余方程）（5）解方程（6）后处理（扩展计算）

3. Matlab有限元分析实战【实例1】

分析：步骤一：单元划分

步骤二：构造单元刚度矩阵 >>k1=SpringElementStiffness(100) >>…?

步骤三：构造系统刚度矩阵 a) 分析SpringAssemble库函数function y = SpringAssemble(K,k,i,j) % This function assembles the element stiffness % matrix k of the spring with nodes i and j into the % global stiffness matrix K. % function returns the global stiffness matrix K % after the element stiffness matrix k is assembled. K(i,i) = K(i,i) + k(1,1); K(i,j) = K(i,j) + k(1,2); K(j,i) = K(j,i) + k(2,1); K(j,j) = K(j,j) + k(2,2); y = K; b) K是多大矩阵？今天的系统刚度矩阵是什么？因为 11 22 1212 k k k k k k k k - ?? ?? - ????--+ ?? 所以 1000100 0200200 100200300 - ?? ?? - ?? ?? -- ?? ？

《有限元基础教程》_【MATLAB算例】3.3.7(2)__三梁平面框架结构的有限元分析(Beam2D2Node)

【MA TLAB 算例】3.3.7(2) 三梁平面框架结构的有限元分析 (Beam2D2Node) 如图3-19所示的框架结构，其顶端受均布力作用，结构中各个截面的参数都为：113.010Pa E =?，746.510I m -=?，426.810A m -=?。试基于MA TLAB 平台求解该结构的节点位移以及支反力。图3-19 框架结构受一均布力作用解答：对该问题进行有限元分析的过程如下。（1）结构的离散化与编号将该结构离散为3个单元，节点位移及单元编号如图3-20所示，有关节点和单元的信息见表3-5。（a ）节点位移及单元编号

（b）等效在节点上的外力图3-20 单元划分、节点位移及节点上的外载（2）各个单元的描述首先在MA TLAB环境下，输入弹性模量E、横截面积A、惯性矩I和长度L，然后针对单元1，单元2和单元3，分别二次调用函数Beam2D2Node_ElementStiffness，就可以得到单元的刚度矩阵k1(6×6)和k2(6×6)，且单元2和单元3的刚度矩阵相同。 >> E=3E11; >> I=6.5E-7; >> A=6.8E-4; >> L1=1.44; >> L2=0.96; >> k1=Beam2D2Node_Stiffness(E,I,A,L1); >> k2=Beam2D2Node_Stiffness(E,I,A,L2); （3）建立整体刚度方程将单元2和单元3的刚度矩阵转换成整体坐标下的形式。由于该结构共有4个节点，则总共的自由度数为12，因此，结构总的刚度矩阵为KK(12×12)，对KK清零，然后两次调用函数Beam2D2Node_Assemble进行刚度矩阵的组装。 >> T=[0,1,0,0,0,0;-1,0,0,0,0,0;0,0,1,0,0,0;0,0,0,0,1,0;0,0,0,-1,0,0;0,0,0,0,0,1] ; >> k3=T'*k2*T; >> KK=zeros(12,12); >> KK=Beam2D2Node_Assemble(KK,k1,1,2);

(完整版)有限元大作业matlab---课程设计例子

有限元大作业程序设计学校：天津大学院系：建筑工程与力学学院专业：01级工程力学姓名：刘秀学号：\\\\\\\\\\\ 指导老师：

连续体平面问题的有限元程序分析 [题目]：如图所示的正方形薄板四周受均匀载荷的作用，该结构在边界上受正向分布压力， m kN p 1=，同时在沿对角线y 轴上受一对集中压力，载荷为2KN ，若取板厚1=t ，泊松比0=v 。 [分析过程]：由于连续平板的对称性，只需要取其在第一象限的四分之一部分参加分析，然后人为作出一些辅助线将平板“分割”成若干部分，再为每个部分选择分析单元。采用将此模型化分为4个全等的直角三角型单元。利用其对称性，四分之一部分的边界约束，载荷可等效如图所示。

[程序原理及实现]：用FORTRAN程序的实现。由节点信息文件NODE.IN和单元信息文件ELEMENT.IN，经过计算分析后输出一个一般性的文件DATA.OUT。模型基本信息由文件为BASIC.IN生成。该程序的特点如下：问题类型：可用于计算弹性力学平面问题和平面应变问题单元类型：采用常应变三角形单元位移模式：用用线性位移模式载荷类型：节点载荷，非节点载荷应先换算为等效节点载荷材料性质：弹性体由单一的均匀材料组成约束方式：为“0”位移固定约束，为保证无刚体位移，弹性体至少应有对三个自由度的独立约束方程求解：针对半带宽刚度方程的Gauss消元法

输入文件：由手工生成节点信息文件NODE.IN，和单元信息文件ELEMENT.IN 结果文件：输出一般的结果文件DATA.OUT 程序的原理如框图：

Matlab有限元分析操作基础共11页

Matlab有限元分析20140226 为了用Matlab进行有限元分析，首先要学会Matlab基本操作，还要学会使用Matlab进行有限元分析的基本操作。 1. 复习：上节课分析了弹簧系统 x 推导了系统刚度矩阵

3. Matlab有限元分析实战【实例1】

分析：步骤一：单元划分

>>k1=SpringElementStiffness(100)

a) 分析SpringAssemble库函数 function y = SpringAssemble(K,k,i,j) % This function assembles the element stiffness % matrix k of the spring with nodes i and j into the % global stiffness matrix K. % function returns the global stiffness matrix K % after the element stiffness matrix k is assembled. K(i,i) = K(i,i) + k(1,1); K(i,j) = K(i,j) + k(1,2); K(j,i) = K(j,i) + k(2,1); K(j,j) = K(j,j) + k(2,2); y = K; b) K是多大矩阵？今天的系统刚度矩阵是什么？因为 11 22 1212 k k k k k k k k - ?? ?? - ????--+ ?? 所以 1000100 0200200 100200300 - ?? ?? - ????-- ??？

基于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)

天然气管道穿孔局部应力应变分析

天然气管道穿孔局部应力应变分析发表时间：2020-03-24T09:49:13.023Z 来源：《文化时代》2020年1期作者：张益 [导读] 本文主要以X70天然气管道为研究对象，针对穿孔管道的局部力学特性进行分析，通过模拟针对穿孔管道的局部等效应力和塑性应变分布状况进行分析。中国石油天然气管道公司中原输油气分公司山东省德州市 253000 摘要：本文主要以X70天然气管道为研究对象，针对穿孔管道的局部力学特性进行分析，通过模拟针对穿孔管道的局部等效应力和塑性应变分布状况进行分析。关键词：天然气管道；穿孔；局部应力；应变引言天然气是一种高效的清洁能源，目前在生产生活中的应用非常广泛，而管道运输是天然气输送的主要方式，这种推广方式具有安全、高效的特征。天然气管道在长期运行过程中不可避免的会受到腐蚀作用影响，腐蚀深度不断增加会最终导致天然气管道出现穿孔现象，进而引发天然气泄漏，造成不可挽回的后果。 1 天然气管道穿孔模型 1.1穿孔实验模型天然气管道在出现腐蚀现象后，随着时间的不断推移，发生腐蚀的位置会逐渐扩散，最终会形成穿孔。本次实验中选择的天然气穿孔管道内壁直径达到20mm，外壁直径为6mm[1]。 1.2穿孔有限元模型以上述天然气穿孔管道模型为基础，充分利用Solid185单元来建立起从内向外以及从外向内两种穿孔管道模型，将管道利用自由网格进行划分，并针对发生穿孔位置附近的管道进行网格加密，并在此基础上对网格质量进行多次性的改善。 1.3材料模型本次研究中主要选取了X70管线钢天然气管道为模型，这种天然气管道材质本身的弹性模量达到了210Gpa，柏松比达到0.3。该管材具备了一定的连续屈服特征，而且没有明显的屈服平台，针对建立模型进行多线推动强化，以此来描述管道本身的弹塑性[2]。 2 内压对最大应力-应变的影响最大 2.1 应力-应变随内压变化分析为了能够针对天然气管道穿孔在不同的压力状况下局部位置的应力以及应变分布状况进行全面分析。针对天然气管道内壁施加了一个压力为25.0MPa的内压，与此同时设置了50个子部，也就是说，每一个子部增加表示内压升高了0.5MPa，针对每一个子部的最终计算结果进行详细统计之后就能够最终得出不同压力状况下天然气管道的应力-应变分布状况。在针对天然气管道穿孔局部最大等效应力、塑性应变变化趋势进行分析，为了能够对其变化状况进行更加清晰的展示，以16.0MPa为基点将所有应变数据划分成两组，并分别绘制曲线。针对最终绘制出的曲线进行分析后可以知道，在最大等效应力、塑性应变变化方面内外穿孔相似度非常高，当内压上升到5.0MPa的情况下，最大应力增长趋势趋于缓慢。而与穿孔位置距离较远的位置开始出现塑性应变时，内压达到了16.0MPa，而此时，天然气管道发生穿孔的位置，最大应力、应变增长速度开始明显变大。之所以出现这种现象是因为只有穿孔位置周围的天然气管道进入了塑性区，其他部分天然气管道管壁仍然处在弹性阶段，而天然气管道的弹性性能对塑性区塑性流动会产生一定的限制作用，导致塑性区实际产生的应变并不明显，而随着整个管道大部分位置进入塑性区之后，穿孔位置附近实际产生的塑性流动受到了限制作用也逐渐减小，在此基础上使得应变出现了明显增加现象。随着内压的进一步增加，达到19.5MPa的时候，穿孔位置的最大应力达到了极限强度，因此开始逐渐趋于稳定。内压进一步增长到20.0MPa的情况下，内外穿孔位置附近最大塑性硬件呈现出指数倍的增长，在这种情况下天然气管道非常容易出现开裂现象。而管道穿孔之后，内压与正常运行压力相比较要小很多，因此要想达到20.0MPa比较困难，因此常温状态下通常不会出现开裂问题。 2.2 应力应变云图分析在针对不同压力条件下穿孔局部应力应变云图技术分析可以知道，在穿孔位置的外壁边缘出现了应力-应变最大值，而且在天然气管道的径向方向上分布着较大的应力-应变。当天然气管道内压达到16.0MPa的情况下，整个天然气管壁开始出现屈服现象，当内压进一步缓慢增加的时候，天然气管道关键部位最大应力应变出现了快速的增加现象，穿孔位置周边较大的应变分布范围也在迅速扩大；当内压达到19.0MPa的情况下，应变值超过0.026的分布范围外边缘与穿孔位置的距离已经非常远；当内压进一步增加，达到20.0MPa的时候，天然气管道的绝大部分管壁的应变值已经超过了0.026，沿着厚度方向天然气管道应变值分布在0.077~0.231和范围内[3]，由此也可以知道，天然气管道的穿孔开裂首先会从关键点开始，对沿着管壁的厚度方向逐渐形成贯穿性裂纹。 3 管道各参数对最大附近应变影响分析 3.1 穿孔尺寸影响当天然气内压在20.0MPa情况下，分析最大应变于穿孔半径的关系趋势可以发现，随着穿孔孔径的逐渐增加，最大应变值在逐渐减小，当穿孔孔径超过一定数值的时候，最大塑性应变波动呈现出复杂化。这主要是因为，当穿孔半径相对比较小的时候，仅仅在穿孔的外壁边缘位置出现最大塑性应变，而当其超过某一个数值时，发生最大塑性应变的位置也会逐渐向着中间移动，这样就导致应变值的变化更加复杂。 3.2 管道壁厚对最大塑性硬件影响天然气管道的壁厚对管道本身承载能力的影响非常大，因此天然气管道穿孔局部应力-应变分布状况也必然会受到管道壁厚的巨大影响。针对内压为20.0MPa情况下不同管道壁厚下最大应变与和穿孔距离较远位置的应变变化趋势分析可以知道。在壁厚不断增加的情况下，穿孔局部最大促进应变会出现明显下降，而且与距离穿孔位置较远位置的管壁应变变化状况相比较，穿孔局部实际发生的最大促进应

matlab有限元分析实例

1.物理现象：这个对工程师来说是直观的物理现象和物理量，温度多少度，载荷是多大等等。通常来说，用户界面中呈现的、用户对工程问题进行设置时输入的都是此类信息。 2.数学方程：将物理现象翻译成相应的数学方程，例如流体对应的是NS方程，传热对应的是传热方程等等；大部分描述这些现象的方程在空间上都是偏微分方程，偶尔也有ODE（如粒子轨迹、化学反应等）。在这个层面，软件把物理现象“翻译” 为以解析式表示的数学模型。 3.数值模型：在定义了数学模型，并执行了网格剖分后，商业软件会将数学模型离散化，利用有限元方法、边界元法、有限差分法、不连续伽辽金法等方法生成数值模型。软件会组装并计算方程组雅可比矩阵，并利用求解器求解方程组。这个层面的计算通常是隐藏在后台的，用户只能通过一些求解器的参数来干预求解。有限元是一种数值求解偏微分方程的方法。基本过程大致是设置形函数，离散，形成求解矩阵，数值解矩阵，后处理之类的。 MATLAB要把这些过程均自己实现，不过在数值求解矩阵时可以调用已有函数。可以理解为MATLAB是一个通用的计算器，当然它的功能远不止如此。

而ANSYS之类的叫做通用有限元软件，针对不同行业已经将上述过程封装，前后处理也比较漂亮，甚至不太了解有限元理论的人也能算些简单的东西，当然结果可靠性又另说了。比较两者，ANSYS之类的用起来容易得多，但灵活性不如MATLAB。MATLAB用起来很困难，也有人做了一些模块，但大多数只能解决一些相对简单的问题。对于大多数工程问题，以及某些领域的物理问题，一般都用通用有限元软件，这些软件还能添加一些函数块，用以解决一些需要额外设置的东西。但是对于非常特殊的问题，以及一般性方程的有限元解，那只能用MATLAB或C，Fortran之类的了。

机械零件的应力应变分析

§3－3机械零件的应力应变分析一、拉（压）杆应力应变分析（一）应力分析前面应用截面法，可以求得任意截面上内力的总和，现在进一步分析横截面上的应力情况，首先研究该截面上的内力分布规律，内力是由于杆受外力后产生变形而引起的，我们首先通过实验观察杆受力后的变形现象，并根据现象做出假设和推论；然后进行理论分析，得出截面上的内力分布规律，最后确定应力的大小和方向。现取一等直杆，拉压变形前在其表面上画垂直于杆轴的直线和（图3-28）。拉伸变形后，发现和仍为直线，且仍垂直于轴线，只是分别平行地移动至和。于是，我们可以作出如下假设：直杆在轴向拉压时横截面仍保持为平面。根据这个“平面假设”可知，杆件在它的任意两个横截面之间的伸长变形是均匀的。又因材料是均匀连续的，所以杆件横截面上的内力是均匀分布的，即在横截面上各点处的正应力都相等。若杆的轴力为，横截面积为,,于是得： ???????????????????????? (3－2) 这就是拉杆横截面上正应力的计算公式。当为压力时，它同样可用于压应力计算。规定拉应力为正，压应力为负。例3－3? 图3-29（a）为一变截面拉压杆件，其受力情况如图示，试确定其危险截面。解? 运用截面法求各段内力，作轴力图[图3-29（b）]：段：????????? 段：段：???????? 段：根据内力计算应力，则得：段：????????? 段：

段：最大应力所在的截面称为危险截面。由计算可知，段和段为危险截面。（二）、拉（压）杆的变形杆件受轴向拉力时，纵向尺寸要伸长，而横向尺寸将缩小；当受轴向压力时，则纵向尺寸要缩短，而横向尺寸将增大。设拉杆原长为，横截面面积为(图3-30)。在轴向拉力P作用下，长度由变为，杆件在轴线方向的伸长为, 。实验表明，工程上使用的大多数材料都有一个弹性阶段，在此阶段范围内，轴向拉压杆件的伸长或缩短量，与轴力和杆长成正比，与横截面积成反比。即，引入比例常数则得到： ??????????????????? （3－3）这就是计算拉伸（或压缩）变形的公式，称为胡克定律。比例常数称为材料的弹性模量，它表征材料抵抗弹性变形的性质，其数值随材料的不同而异。几种常用材料的值已列入表3－1中。从公式（3－3）可以看出，乘积越大，杆件的拉伸（或压缩）变形越小，所以称为杆件的抗拉（压）刚度。上式改写为：其中，而表示杆件单位长度的伸长或缩短，称为线应变(简称应变)，即。是一个无量纲的量，规定伸长为正，缩短为负。则（3－3）式可改写为：????????????????????????????????????????????? ?????????????????????????????????????????????????????? （3-4）式（3-17）表示，在弹性范围内，正应力与线应变成正比。这一关系通常称为单向胡克定律。杆件在拉伸（或压缩）时，横向也有变形。设拉杆原来的横向尺寸为，变形后为（图3-30），则横向应变为：实验指出，当应力不超过比例极限时，横向应变与轴向应变之比的绝对值是一个常数。即称为横向变形系数或泊松比，是一个无量纲的量。和弹性模量E一样，泊松比也是材料固有的弹性常数。因为当杆件轴向伸长时，横向缩小；而轴向缩短时，横向增大，所以和符号是相反的。

Matlab-PDE工具箱有限元法求解偏微分方程

在科学技术各领域中，有很多问题都可以归结为偏微分方程问题。在物理专业的力学、热学、电学、光学、近代物理课程中都可遇见偏微分方程。偏微分方程，再加上边界条件、初始条件构成的数学模型，只有在很特殊情况下才可求得解析解。随着计算机技术的发展，采用数值计算方法，可以得到其数值解。偏微分方程基本形式而以上的偏微分方程都能利用PDE工具箱求解。 PDE工具箱 PDE工具箱的使用步骤体现了有限元法求解问题的基本思路，包括如下基本步骤： 1) 建立几何模型 2) 定义边界条件 3) 定义PDE类型和PDE系数 4) 三角形网格划分

5) 有限元求解 6) 解的图形表达以上步骤充分体现在PDE工具箱的菜单栏和工具栏顺序上，如下具体实现如下。打开工具箱输入pdetool可以打开偏微分方程求解工具箱，如下首先需要选择应用模式，工具箱根据实际问题的不同提供了很多应用模式，用户可以基于适

当的模式进行建模和分析。在Options菜单的Application菜单项下可以做选择，如下或者直接在工具栏上选择，如下列表框中各应用模式的意义为： ① Generic Scalar：一般标量模式（为默认选项）。 ② Generic System：一般系统模式。 ③ Structural Mech.，Plane Stress：结构力学平面应力。 ④ Structural Mech.，Plane Strain：结构力学平面应变。

⑤ Electrostatics：静电学。 ⑥ Magnetostatics：电磁学。 ⑦ Ac Power Electromagnetics：交流电电磁学。 ⑧ Conductive Media DC：直流导电介质。 ⑨ Heat Tranfer：热传导。 ⑩ Diffusion：扩散。可以根据自己的具体问题做相应的选择，这里要求解偏微分方程，故使用默认值。此外，对于其他具体的工程应用模式，此工具箱已经发展到了Comsol Multiphysics软件，它提供了更强大的建模、求解功能。另外，可以在菜单Options下做一些全局的设置，如下 l Grid：显示网格 l Grid Spacing…：控制网格的显示位置 l Snap：建模时捕捉网格节点，建模时可以打开 l Axes Limits…：设置坐标系范围 l Axes Equal：同Matlab的命令axes equal命令建立几何模型使用菜单Draw的命令或使用工具箱命令可以实现简单几何模型的建立，如下各项代表的意义分别为

弹性力学基础分析

岩石力学－第三讲：弹性力学基础（一、应力应变分析）教学备忘录

大多数物质在受到外力时发生变形，在外力撤除后又能恢复到原来的形状。我们把物质的这种性质称之为弹性。弹性是岩石力学的基础，外力和相应的变形间呈线性关系是最简单的情况。当在外力的作用下，物质发生的变形足够小，那么这种关系几乎总是线性的。因此，线弹性是所有弹性问题的基础。1.1介绍了固体物质的线弹性特性。在实际情况下，线弹性的有效区域经常被超越。1.1中介绍了一些岩石非线性行为的一般特征。在石油工程岩石力学中，更多的兴趣集中在那些具有有效孔隙和渗透性的岩石上。固体材料的弹性理论不能完全描述这种介质，因此，应该引入多孔弹性的概念。岩石的弹性反应也可能是与时间相关的，因此，介质的变形也是随着时间而变化的，甚至在外力不变的情况下也是这样。1.3节和1.4节分别介绍了多孔物质的弹性特性和随时间变化效应。 1.1 线性弹性理论弹性理论建立在应力和应变这两个概念之上，在1.1.1和1.1.2节中对应力和应变分别做了介绍。1.1.3节和1.1.4节分别介绍了各向同性介质和各向异性介质应力和应变之间的线性本构方程 1.1.1 应力考虑图1.1所示（见多媒体）的情况，一个重物加在柱子的顶部。由于重物的重量，一个作用力施加在柱子上，同时柱子会产生一个大小相等、方向相反的力。而柱子本身支撑在地面上，因此，施加在柱子顶部的作用力必然会通过柱子的任意横截面。 a)处横截面的区域如A 所示。如果施加在横截面上的力为F ，则该截面处的应力σ定义为： A F = σ （1.1）应力经常用Pa(=Pascal=N/m 2 )、 bar 、atmosphere 、 psi(=lb/sq.inch.)或 dynes/cm 2 等单位来表示。在理论计算中，国际单位Pa 是最合适的单位，而其它单位大多应用于工程计算。应力符号σ不仅表示受力面的物理性质，而且已经依照惯例进行了定义。在岩石力学中，符号惯例规定：压应力为正。历史原因在于：岩石力学涉及到的应力几乎都是压应力。当符号惯例被一直使用时并没有引发问题，但是，记住一些其它科学，包括弹性力学使用相反的符号惯例是重要的。正如公式（1.1）所表明的那样，应力被一个力和一个截面（或通常来说是一个平面）所定义，力是被施加的。看看b)处的截面，施加在截面上的力等于施加在截面a)处的力（忽略柱子本身的重量）。然而，b)处横截面的区域A ′明显小于A 。因此，b)处的应力σ′＝F/A ′大于a)处的应力，即在受力试件中，应力随位置变化而变化。我们可以将a)处截面分为无数个小单元ΔA ，总力F 的一个无限小单元力ΔF 施加在这个小单元ΔA 上（图1.2）。不同的小单元，力ΔF 也不同。设想一小单元i ，其包含一点P 。当其面积Δai 趋近于零时，点P 处的应力被定义为Δfi/Δai 的极限，即： i i A A F i ??=→?lim σ （1.2）

有限元钢架结构分析手算matlabansys模拟

有限元大作业——钢架结构分析选题人：日期：2016年6月2日

目录：第一章：问题重述 (2) 一、题目内容： (3) 二、题目要求： (3) 第二章：有限元法手工求解 (3) 一、平面两单元离散化 (4) 二、单元分析 (4) 三、单元组装 (6) 四、边界条件引入及组装总体方程 (7) 五、求解整体刚度方程，计算节点2的位移和转角 (7) 六、求节点1、3支撑反力 (8) 七、设定数据，求解结果 (8) 八、绘制轴力图、弯矩图、剪力图 (9) 第三章、matlab编程求解： (11) 一、总体流程图绘制： (11) 二、输入数据： (12) 三、计算单元刚度矩阵： (12) 四、建立总体刚度矩阵： (13) 五、计算未约束点位移： (13) 六、计算支反力： (13) 七、输出数据： (13) 八、编程： (13) 第四章有限元求解 (13) 一、预处理 (13) 二、模型建立： (15) 二、分析计算 (17) 三、求解结果 (18) 四、绘制图像 (19) 第五章结果比较 (22) 第六章心得体会 (22) 第七章附录 (23) 一、matlab程序 (24) 第一章：问题重述

一、题目内容：图示平面钢架结构图题目内容二、题目要求：（1）采用平面梁单元进行有限元法手工求解，要求写出完整的求解步骤，包括： a）离散化：单元编号、节点编号； b）单元分析：单元刚度矩阵，单元节点等效载荷向量； c）单元组长：总体刚度矩阵，总体位移向量，总体节点等效载荷； d）边界条件的引入及总体刚度方程的求解； e）B点的位移，A、C处支撑反力，并绘制该结构的弯矩图、剪力图和轴力图。（2）编制通用平面钢架分析有限元Matlab程序，并计算盖提，与手工结果进行比较；（3）利用Ansys求解，表格列出B点的位移，A、C处支反力，绘制弯矩图、剪力图和轴力图，并与手算和Matlab程序计算结果比较。（4）攥写报告，利用A4纸打印；（5）心得体会，并简要说明各成员主要负责完成的工作。第二章：有限元法手工求解

利用Matlab进行有限元分析结果的可视化显示

利用Matlab进行有限单元法计算结果的可视化显示摘要本文用一个简单的例子给出了用Matlab进行有限单元法计算结果可视化显示的方法。采用Matlab进行可视化显示，可以在获得较好的可视化显示效果的基础上，节省科研人员的大量时间和精力。关键字：有限元，后处理，可视化，Matlab 有限单元法是工程数值分析的有力工具，可以应用于固体力学、结构分析、温度场模拟等诸多领域。有限单元法一般可以分为前处理、计算以及后处理三部分，市场上现有的有限元商业软件都提供了这三部分功能模块。但有时，由于各种原因，科研人员必须自行编写有限元分析程序，作者通过自身实践，认为Matlab可以较好的进行有限单元法计算结果的可视化显示。 Matlab由美国MathWorks公司开发，历经二十多年的发展，现已成为国际公认的优秀科技应用软件之一，在机械、航天、医药等多个科研、工程领域有着广泛的应用。Matlab 本身具有丰富的可视化显示手段，但遗憾的是，目前对于Matlab的应用研究主要集中在其强大的科学计算能力方面，而对科学计算结果的可视化显示，尤其对由空间点云构成的形体的可视化显示研究涉及甚少，作者通过查阅相关资料，以及探索和实践，成功地进行了三维形体有限元分析结果的可视化显示。 1.准备数据针对Matlab对空间点云构成形体的数据格式要求，必须重新编排有限元分析中前处理部分以及计算部分所获得的数据。下面以空间单位立方体为例，介绍Matlab对数据文件格式的要求。若有空间单位正方体，将其划分为四面体网格，图1为该正方体的节点编号及其网格拓朴结构，表1为节点的坐标值以及节点处的有限元计算结果（此处为温度）。表1：单位正方体顶点坐标及其温度图1：空间立方体顶点编号及其网格拓朴结构

有限元的MATLAB解法

有限元的MATLAB解法 1.打开MATLAB。 2.输入“pdetool”再回车，会跳出PDE Toolbox的窗口（PDE意为偏微分方程，是partial differential equations的缩写），需要的话可点击Options菜单下Grid命令，打开栅格。 3.完成平面几何模型：在PDE Toolbox的窗口中，点击工具栏下的矩形几何模型进行制作模型，可画矩形R，椭圆E，圆C，然后在Set formula栏进行编辑并（如双脊波导R1+R2+R3改为RI-R2-R3,设定a、b、s/a、d/b的值从而方便下步设定坐标）用算术运算符将图形对象名称连接起来，若还需要，可进行储存，形成M文件。 4.用左键双击矩形进行坐标设置：将大的矩形left和bottom都设为0，width是矩形波导的X轴的长度，height是矩形波导的y轴的长度，以大的矩形左下角点为原点坐标为参考设置其他矩形坐标。 5.进行边界设置：点击“Boundary”中的“Boundary Mode”,再点击

“Boundary”中的“Specify Boundary Conditions”,选择符合的边界条件，Neumann为诺曼条件,Dirichlet为狄利克雷条件，边界颜色显示为红色。 6.进入PDE模式：点击"PDE"菜单下“PDE Mode”命令，进入PDE 模式，单击“PDE Specification”,设置方程类型，“Elliptic”为椭圆型，“Parabolic”为抛物型，“Hyperbolic”为双曲型，“Eigenmodes”为特征值问题。 7.对模型进行剖分：点击“Mesh”中“Initialize Mesh”进行初次剖分，若要剖的更细，再点击“Refine Mesh”进行网格加密。 8.进行计算：点击“Solve”中“Solve PDE”,解偏微分方程并显示图形解，u值即为Hz或者Ez。 9.单击“Plot”菜单下“Parameters”选项，打开“Plot Selection”对话框。选中Color,Height(3-D plot)和Show mesh三项，然后单击“Plot”按钮，显示三维图形解。 10.如果要画等值线图和矢量场图，单击“Plot”菜单下“Parameters”选项,打开“Plot Selection”对话框。选中Contour和Arrows两项，然后单击Plot按钮，可显示解的等值线图和矢量场图。 11.将计算结果条件和边界导入MATLAB中：点击“Export Solution”,再点击“Mesh”中“Export Mesh”。

局部应力应变法

局部应力应变法传统的局部应力应变法以Manson 一Coffin 公式为材料疲劳性能曲线.以应力集中处的局部点应力作为衡量结构受载严重程度的参数.这一方法在大应变低寿命时与实际情况符合很好.但进人高周疲劳，由于Manson 一Coffin 公式与实验结果的差距逐渐增大，由于缺口根部塑性的消失而使应力梯度变大，致使传统的局部应力应变法过低地估计了结构的疲劳寿命.就实际工程结构而育，通常受到随机载荷的作用，在大多数情况下，载荷谱中的高载处于低周疲劳阶段，大多数的中低级载荷处于高周疲劳阶段，所以寻找一个同时适用于高周和低周疲劳寿命估算的方法是其有很大实际意义的。 ( ε-f N ) 曲线是是重要的材料疲劳性能曲线，在局部应力应变法中，它是结构疲劳寿命估算的基本性能数据。传统的局部应力应变法采用Manson-Coffin 公式来描述 ''(2)(2)f b c a f f f N N E σεε=+ (1) Manson-Coffin 公式虽然在工程上得到了广泛的应用，但也存在着一些严重的不足：①大多数金属材料按Manson-Coffin 分解后的塑性线不能很好地用直线来拟合，而是向下弯曲的曲线；②Manson-Coffin 公式仅适用于解决低周疲劳寿命的计算，而在高周疲劳时计算出的寿命与实验结果相差较大；③当（1）式中的 f N 趋于无穷时，ε趋于零，即Manson-Coffin 公式没有反映出的疲劳极限，这与实际情况不符。文献[1]针对传统的局部应力应变法存在的这两个缺陷，提出解决这一问题的方法:用等效应变一寿命曲线或四参数应变一寿命曲线替换Manson 一Coffin 公式，用更合适的缺口疲劳系数或缺口场强度来描述缺口受载的严重程度，希望将传统的局部应力应变法推广到高周疲劳寿命的估算。四参数（ε-f N ）曲线：在中高疲劳区（1）式已不太适用，文献[2]提出了一个四参数的（ε-f N ）曲线拟合公式 2013lg(/)lg *ln{}lg(/) t f t A N A A A εε?=+? (2) 式中：为四个回归参数。（2）式具有以下特点：①它适用于中高周疲劳阶段（ε-f N ），克服了Manson-Coffin 公式在高周疲劳段误差较大的缺点；②当 f N 趋于无穷时，ε趋近于A3，可反映出材料的疲劳极限；③与当量应变-寿命曲线公式（1）相比，不需由实验给出参数同时比（ε- f N ）曲线的拟合度高。文献[3]指出现行方法在计算中低周疲劳有较好的寿命预测精度，但对高周疲劳寿命预测精度不高。它认为主要是因为没有考虑到应力集中、表面加工状况、尺寸和环境介质的

有限单元法matlab编程实例

主程序 E=210e6;A=2e-2;I=5e-5;L1=3;L2=4;L3=3; k1=PlaneFrameElementStiffness(E,A,I,L1,90); k2=PlaneFrameElementStiffness(E,A,I,L2,0); k3=PlaneFrameElementStiffness(E,A,I,L3,270); K=zeros(12,12); K=PlaneFrameAssemble(K,k1,1,2); K=PlaneFrameAssemble(K,k2,2,3); K=PlaneFrameAssemble(K,k3,3,4) k=K(4:9,4:9); f=[-20;0;0;0;0;12]; u=k\f U=[0;0;0;u;0;0;0] F=K*U u1=[U(1);U(2);U(3);U(4);U(5);U(6)]; u2=[U(4);U(5);U(6);U(7);U(8);U(9)]; u3=[U(7);U(8);U(9);U(10);U(11);U(12)]; f1=PlaneFrameElementForces(E,A,I,L1,90,u1) f2=PlaneFrameElementForces(E,A,I,L2,0,u2) f3=PlaneFrameElementForces(E,A,I,L3,270,u3)

需调用的函数和子程序 function y=PlaneFrameAssemble(K,k,i,j) %PlaneFrameAssemble This function assembles the element stiffness %matrix k of the plane frame element with nodes i and j into the global %stiffness matrix K .This function returns the global stiffness matrix K after %the element stiffness matrix k is assembled. K(3*i-2,3*i-2)=K(3*i-2,3*i-2)+k(1,1); K(3*i-2,3*i-1)=K(3*i-2,3*i-1)+k(1,2); K(3*i-2,3*i)=K(3*i-2,3*i)+k(1,3); K(3*i-2,3*j-2)=K(3*i-2,3*j-2)+k(1,4); K(3*i-2,3*j-1)=K(3*i-2,3*j-1)+k(1,5); K(3*i-2,3*j)=K(3*i-2,3*j)+k(1,6); K(3*i-1,3*i-2)=K(3*i-1,3*i-2)+k(2,1); K(3*i-1,3*i-1)=K(3*i-1,3*i-1)+k(2,2); K(3*i-1,3*i)=K(3*i-1,3*i)+k(2,3); K(3*i-1,3*j-2)=K(3*i-1,3*j-2)+k(2,4); K(3*i-1,3*j-1)=K(3*i-1,3*j-1)+k(2,5); K(3*i-1,3*j)=K(3*i-1,3*j)+k(2,6); K(3*i,3*i-2)=K(3*i,3*i-2)+k(3,1); K(3*i,3*i-1)=K(3*i,3*i-1)+k(3,2); K(3*i,3*i)=K(3*i,3*i)+k(3,3); K(3*i,3*j-2)=K(3*i,3*j-2)+k(3,4); K(3*i,3*j-1)=K(3*i,3*j-1)+k(3,5); K(3*i,3*j)=K(3*i,3*j)+k(3,6); K(3*j-2,3*i-2)=K(3*j-2,3*i-2)+k(4,1); K(3*j-2,3*i-1)=K(3*j-2,3*i-1)+k(4,2); K(3*j-2,3*i)=K(3*j-2,3*i)+k(4,3); K(3*j-2,3*j-2)=K(3*j-2,3*j-2)+k(4,4); K(3*j-2,3*j-1)=K(3*j-2,3*j-1)+k(4,5); K(3*j-2,3*j)=K(3*j-2,3*j)+k(4,6); K(3*j-1,3*i-2)=K(3*j-1,3*i-2)+k(5,1); K(3*j-1,3*i-1)=K(3*j-1,3*i-1)+k(5,2); K(3*j-1,3*i)=K(3*j-1,3*i)+k(5,3); K(3*j-1,3*j-2)=K(3*j-1,3*j-2)+k(5,4); K(3*j-1,3*j-1)=K(3*j-1,3*j-1)+k(5,5); K(3*j-1,3*j)=K(3*j-1,3*j)+k(5,6); K(3*j,3*i-2)=K(3*j,3*i-2)+k(6,1); K(3*j,3*i-1)=K(3*j,3*i-1)+k(6,2); K(3*j,3*i)=K(3*j,3*i)+k(6,3); K(3*j,3*j-2)=K(3*j,3*j-2)+k(6,4); K(3*j,3*j-1)=K(3*j,3*j-1)+k(6,5); K(3*j,3*j)=K(3*j,3*j)+k(6,6); y=K;