MATLAB数值模拟.ppt

1、MATLAB在数值模拟中的应用,报告人：海洋孤树,提纲,1.Matlab一些常见有限元开源 程序包的简单介绍 2.Matlab与其他软件（包）的结合 3.光滑有限元,1.Matlab常见有限元程序包,Matalb是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际

2、科学计算软件的先进水平。 Matlab的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用Matlab来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多。因而Matlab用来验证一种新的数值方法是十分方便的。,1.1 Matlab简介,2.1 Matlab编写的有限元程序及程序包,1） 最经典的是Jack Chessa编写的有限元程序包，详细介绍见Programing the Finite Element Method with Matlab。这个程序包主要是常规有限元2维弹性问题的一些程序，它包括前后处理程序及常见的如四节点、八节点、九节点等参单元，三节点三

3、角形单元和六节点三角形单元，是有限元入门学习的工具。,2) MATLAB Codes for Finite Element Analysis 这是Ferreira, A. J. M编写的Matlab关于有限元的书，其中包括2维、3维桁架机构，框架结构，Mindlin 板及层合板的分析，模块性好，对研究板壳问题初学者提供了很好的指导作用。,3). AFEMMATLAB,这是一个能够进行网格自适应的有限元程序包，详细参照AFEMMATLAB: A MATLAB PACKAGE OF ADAPTIVE FINITE ELEMENT METHODS,4). mlfem,这是一个运用于生物力学分析的程序

4、包。主要用于模拟骨骼损伤分析。,5). Openfem 这是一个结构化的有限元程序包。能够和c、fortran语言结合编程提高计算效率，另外该程序包可以进行各向异性材料，非线性材料、板壳问题的分析，包括2维，3维问题。该程序包还包括了多场耦合单元，必要压电单元。,6）.FAESOR 该程序包是由Petr Krysl课题组编写的Matlab面对对象的有限元程序包，该程序包一直都在更新，最新版本更新到了2012年4月13日 。 该程序包采用面向编程方法，程序效率较高，本身带有生成复杂网格的子程序包。应用范围主要包括接触分析，不可压缩材料分析、电热分析，热分析、声学分析、波动分析、弹塑性分析、超弹性

5、材料分析、动力分析等等。,2.Matlab与其他软件包的结合,2.1下面介绍Matlab与Ansys结合 Matlab编程简便，在现有程序包基础上能够应用于分析新的问题，新方法，比Ansys二次开发灵活。然而matlab本身不具备很好的网格划分功能，这需要调用,首先在Ansys中进行建模、网格划分，编写ansys.mac文件读取网格单元，节点，边界信息.然后在matlab中编写网格读取文件读取网格信息。 ANSYS2Matlab是一个开源的工具包，该工具包括了Ansys.mac、matlab读取ansys网格信息的文件。主要步骤是： 1.添加ansys_ch.mac、start120.ans至

6、ansys路径，打开Ansys软件，在ansys工具栏将会出现ansys_ch按钮用于保存ansys网格节点，单元，边界信息。,2. 通过loadFromGridFile.m读取Ansys网格信息，并在matlab中重新生成网格图。,右图是一个简单带圆孔平面板问题，取1/4并网格划分。用Ansys2Matlab工具使得Matlab成功调用Ansys网格文件,Ansys网格,Matlab重生成的网格,Ansys2Matlab工具仅仅能够调用三角形三节点单元、四节四边形单元，且限于2维问题，需要对进行修改一边应用高阶单元及3为维问题，下面在其基础上修改应用到3维的例子。,Matlab重新生成的网格

7、,Ansys网格,类似地，我们可以通过编写相关程序调用Abaqus、Hypermesh等成熟商业软件的网格文件。,2.2 Matlab调用Lapack程序包,我们知道Lapack是一个非常经典的线性代数程序包，由Fortran编程语言写就。而Matlab通过书写不同的Mex文件可以调用C/C+、Fortran其他编程语言的程序。这样可以使得Matlab计算速度、精度提高。 在编译mex文件前需要进行mex setup操作。在不同操作系统中，我们可能需要对相应的mex文件进行修改。下面举例说明调用Lapack 必要性。 在进行非常怪状的广义矩阵特征值求解时，直接调用Matlab中eig函数，可能

8、出现inf的特征值，这往往与真实解不符。本人通过调用Lapack求解却能够避免这样的问题。,矩阵A：,矩阵B：,下面分别直接应用eig函数和调用lapack中函数zggevx求的广义特征值，可以看出调用lapack中zeggevx求得的特征值和直接使用lapack求得特征值十分接近，而直接使用matlab中eig函数求得特征值相差甚远。,Eig函数求得的特征值：,调用zggevx函数求得的特征值：,Lapack直接求得的特征值：,3.光滑有限元,光滑有限元是刘桂荣教授课题基于应变光滑技术提出来一种新的有限元方法。该方法利用应变光滑对进度矩阵进行改进，无需计算形函数倒数，只需对单元边界进行积分，

9、具有很好的收敛性，对网格畸变影响小等优点。国内湖南大学李光耀课题组已将其应用于汽车碰撞分析。光滑有限元现阶段主要分为基于单元内分块光滑有限元（CSFEM）、边界光滑有限元（ESFEM）、节点光滑有限元、alpha光滑有限元（aSFEM）、面光滑有限元（FSFEM）几大类。,CSFEM示意图,NSFEM示意图,ESFEM示意图,FSFEM示意图,aSFEM示意图,分块光滑有限元原理：,如图所示，二维弹性问题平 衡方程及边界条件表达如下：,为整个问题区域,为边界，其,满足,。,为应力分量，,体力分量。,为界,上牵引力,为界,上位移，,表示为移梯度，,为材料弹张量。,二维弹性体示意图,变分弱形式如下

10、表示：,如图所示，在应变光滑中，设：,为了计算方便，这里取一个简单函数：,由分部积分得：,其中,为一光滑函数，,这里，,为光滑区域。,，,，,由上面公式可以推导出：,类似地，可以推导出光滑应变的公式：,其中光滑应变矩阵,表达如下：,其中，,为光滑区域，,为光滑区域边界，,取高斯点计算只需去光滑区域边界,中点积分即可，则,说明：对于动力问题，质量矩阵求解和常规有限元一样求得,我们在计算单元刚度矩阵时， 使用如图所示积分方式即可。,其他光滑有限元可以采用类似的 方法推导得到。,单元节点形函数值分别为 (1,0,0,0),(0,1,0,0) (0,0,1,0),(0,0,0,1) 其他积分点形函数值可以通过 简单的线性插值得到。,3d CSFEM,算例：,由于光滑有限元基于光滑区域边界积分，对网格畸变影响较小，,甚至对于凹变形单元仍然有效。下面是一个悬臂梁算例,这是三种网格下y=0所有节点