文档库 最新最全的文档下载
当前位置:文档库 › 基于MATLAB的滑坡三维有限元前处理研究

基于MATLAB的滑坡三维有限元前处理研究

基于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有限元分析与应用》由简单到复杂,循序渐进地介绍了各种有限元及其分析与应用方法。书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题,具有很高的工程应用价值。

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 推导了系统刚度矩阵

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的有限元法分析平面应力应变问题刘刚

姓名:刘刚学号: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)

matlab有限元分析实例

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

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

PDE的Galerkin和有限元的MATLAB程序

%Galerkin方法 clear all clc syms x; A=zeros(4,4); for i=1:4 for j=1:4 phy1=x^i; phy2=x^j; dphy1=diff(phy1,1); dphy2=diff(phy2,1); phy=pi^2*phy1*phy2; dphy=dphy1*dphy2; A(i,j)=int(phy+dphy,x,0,1); end end D=[]; for k=1:4 f1=2*pi^2*sin(pi*x)+pi^3*x; f2=x^k; f=f1*f2; D(k)=int(f,x,0,1); end D=D'; C=A\D; C=C'; X=linspace(0,1,6); F=0; for i=1:4 F=F+C(i)*x^i; end for j=1:6 Y(j)=subs(F,X(j)); end Y=Y-pi.*X; Y1=sin(pi.*X); err=norm(abs(Y-Y1)); disp('数值解') disp(Y) disp('整体误差') disp(err)

%%%%%%%%%追赶法 function x=chase(a,b,c,d) n=length(b); u(1)=c(1)/b(1); q(1)=d(1)/b(1); for i=2:n-1 h(i)=b(i)-u(i-1)*a(i-1); u(i)=c(i)/h(i); q(i)=(d(i)-q(i-1)*a(i-1))/h(i); end q(n)=(d(n)-q(n-1)*a(n-1))/(b(n)-u(n-1)*a(n-1)); x(n)=q(n); for i=n-1:-1:1 x(i)=q(i)-u(i)*x(i+1); end end %有限元法 clear all clc m=5; xspan=[0 1]; a=xspan(1);b=xspan(2);h=(b-a)/m; x=[a:h:b]; f=zeros(m-1,1); F1=@(t)(-1/h+(pi^2)*t.*(1-t)*h); c1=quadl(F1,0,1); F2=@(t)(2/h+(pi^2)*(t.^2+(1-t).^2)*h); c2=quadl(F2,0,1); A=diag(c2*ones(m-1,1))+diag(c1*ones(m-2,1),1)+... diag(c1*ones(m-2,1),-1); for i=2:m F=@(t)h*t.*(2*pi^2*sin(pi*(x(i-1)+h*t)))+... h*(1-t).*(2*pi^2*sin(pi*(x(i)+h*t))); f(i-1)=quadl(F,0,1); end v=chase(c1*ones(m-2,1),c2*ones(m-1,1),c1*ones(m-2,1),f); u=[0,v,0]

有限元的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”。

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的命令或使用工具箱命令可以实现简单几何模型的建立,如下 各项代表的意义分别为

matlab有限元解二维抛物方程

%%%%% 真解u=sin(pi*x)*sin(pi*y)*sin(t) %%%%% 方程diff(u,t)-Laplace(u)=f %%%%% f=sin(pi*x)*sin(pi*y)*cos(t)+2*pi^2*sin(pi*x)*sin(pi*y)*sin(t) %clear all % clc %%%%finite element code for parabolic equation with constant coefficient %%%mesh%% node=[0,0;1,0;1,1;0,1]; elem=[2,3,1;4,1,3]; T=1; bdEdge=setboundary(node,elem,’Dirichlet’); n=input(‘Please input initial mesh:’); M=input(‘M=’); for i=1:n [node,elem,bdEdge]=uniformrefine(node,elem,bdEdge); end N=size(node,1); NT=size(elem,1); S=1/NT; r=1/M; A=zeros(N,N); u=zeros(N,M+1); u1=zeros(N,1); f=inline(‘sin(pi*xx(1,1))*sin(pi*xx(1,2))*cos(t)+2*pi^2*sin(pi*xx(1,1))*sin(pi*xx(1,2))*sin(t)’,’xx’,’t’); [lambda,weight]=quadpts(5); p=node’; T=elem’; totalEdge=[elem(:,[2,3]);elem(:,[3,1]);elem(:,[1,2])]; isBdEdge=reshape(bdEdge,3*NT,1); Dirichlet=totalEdge(isBdEdge==1),:); isBdNode=false(N,1); isBdNode(Dirichlet)=true; bdNode=find(isBdNode); freeNode=find(~isBdNode); for j=2:M+1 for i=1:NT F=zeros(N,1); F_ele=zeros(1,3); T_ele=elem(i,:); for m=1:7 xx(m,1)=(p(1,t(1,i))-p(1,t(3,i)))*lambda(m,1)+(p(1,t(2,i))-p(1,t(3,i))*lambda(m,2)+p(1,t(3,i)); xx(m,1)=(p(2,t(1,i))-p(2,t(3,i)))*lambda(m,1)+(p(2,t(2,i))-p(2,t(3,i))*lambda(m,2)+p(2,t(3,i)); end

有限元钢架结构分析手算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编程实例

主程序 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;

Matlab在有限元中的应用

基于MATLAB的有限元结构分析 王剑 (重庆交通大学土木建筑学院,重庆400074) 摘要: MATLAB 是当今国际科学界最具影响力和活力的软件。文章介绍了MATLAB 语言的特点,详细介绍了用MATLAB语言编写结构内力的有限元方法,并通过实例对平面钢架结构进行了内力分析。 关键词:MATLAB 有限元结构 0引言 MATLAB是Mathworks公司推出的,集算法开发、数值运算、符号运算以及图形处理等强大功能于一体的高级技术计算语言和交互式环境。MATLAB意为矩阵实验室,它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中,为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案,并在很大程度上摆脱了传统非交互式程序设计语言的编辑模式,代表了当今国际科学计算软件的先进水平。 有限元法的基本思想是将物体(即连续求解域)离散成有限个且按一定方式相互连接在一起的单元组合,来模拟和逼近原来的物体,从而将一个连续的无限自由度问题简化为离散的有限自由度问题求解的数值分析法。有限元法还有一个特点是,它的理论采用矩阵形式表达。这并不利于一般的计算机语言编制计算机程序,因为传统的计算机语言处理的对象是标量,使用矩阵形式的有限元理论时,必须把矩阵形式的公式转换成标量表示的公式。而如果采用MATLAB,这个特点就变成了有限元法的优点,运算更便捷。 1运用MATLAB编写有限元程序的操作步骤 1.1建立有限元模型 建立有限元模型就是为求解有限元模型做铺垫。需要对节点、单元以及材料的定义。同时对约束条件、集中力、分布力进行定义。然后在M函数文件中以矩阵或向量的形式输入单元号、节点数、材料的性质、约束条件、集中力、分布力。 1.2求解有限元模型 用MATLAB写出每个单元的单元刚度矩阵。按照刚度集成法,把各个单元的刚度矩阵分别放到整体刚度矩阵中的相应位置上,然后根据边界条件进行修正得到整体刚度矩阵。将非节点荷载按照虚功原理分配到单元节点上,从而计算出等效节点力。处理约束条件,最后求解出方程组。 1.3显示计算结果 运行程序后得到各个单元两端的位移列阵和节点力列阵,包括各个单元两端的X方向上的位移、Y方向上的位移和转角,得到的节点力为轴力、剪力、和弯矩。 用MATLAB求解有限元结构内力的算法框图如图1所示。

有限元ANSYS MATLAB 应用

有限元在ANSYS和MATLAB中的应用 工程学院 摘要: 文章简述了有限元分析的基本理论及其求解问题的基本步骤, 介绍了ANSYS 软件的应用,介绍了Matlab 语言特点,给出了Matlab 环境下实现有限元的步骤。说明如何使用Matlab 进行有限元分析,使用该方法进行分析具有精度高、简便、快速及可视化等诸多优点,具有较强的使用价值。 关键词: 有限元分析; ANSYS 软件; 用Matlab 进行有限元分析的优点 1 有限元分析基本理论 有限元分析的基本概念是用较简单的问题代替复杂问题后再求解。它将求解域看成是由许多称为有限元的小的互连子域组成, 对每一单元假定一个合适的近似解, 然后推导求解这个域的满足条件, 从而得到问题的解。这个解不是准确解,而是近似解, 因为实际问题被较简单的问题所代替。由于大多数实际问题难以得到准确解, 而有限元不仅计算精度高, 而且能适应各种复杂形状, 因而成为行之有效的工程分析手段。有限元是那些集合在一起能够表示实际连续域的离散单元[1]。有限元的概念早在几个世纪前就已产生并得到了应用, 例如用多边形逼近圆来求得圆的周长, 但 作为一种方法而被提出, 则是最近的事。有限元法最初被称为矩阵近似方法, 应用于航空器的结构强度计算, 并由于其方便性、实用性和有效性而引起从事力学研究的科学家的浓厚兴趣。经过短短数十年的努力, 随着计算机技术的快速发展和普及, 有限元方法迅速从结构工程强度分析计算扩展到几乎所有的科学技术领域, 成为一种丰富多彩、应用广泛并且实用高效的数值分析方法。有限元方法与其他求解边值问题近似方法的根本区别在于它的近似性仅限于相对小的子域中。 1.1有限元求解问题的分析过程 第一步, 问题及求解域定义: 根据实际问题近似确定求解域的物理性质和几何区域。 第二步, 求解域离散化: 将求解域近似为具有不同有限大小和形状且彼此相连的有限个单元组成的离散域, 习惯上称为有限元网络划分。显然单元越小( 网络越细) 则离散域的近似程度越好, 计算结果也越精确, 但计算量及误差都将增大, 因此求 解域的离散化是有限元法的核心技术之一。 第三步, 确定状态变量及控制方法: 一个具体的物理问题通常可以用一组包含问题状态变量边界条件的微分方程式表示, 为适合有限元求解, 通常将微分方程化为等

matlab有限元分析实例

matlab有限元分析实例 问题描述: 考虑一平面有界区域,设其边界为[。我们求解泊松方程之狄利克雷边值问题。问题的强形式为一椭圆型偏微分方程 当之几何形状稍复杂时,一般无法求得其解析解。我们可应用有限元法来求其数值解。通常我们先写出该问题的抽象弱形式:求使得 其中为检验函数为一适当索伯列夫空间(在本例中也是一希尔伯特空间)为一双线性型为一线性型。其具体表达式为 有限元空间离散 我们采用最简单的二维单元离散单元即三节点线性三角形单元,其插值基函数(即形函数为一次多项式。有限元之核心思想为:使用离散的函数空间来分片逼近连续的函数空间。于是所求之近似解可写成基函数之线性组合 其中为待求系数,常称为自由度。将此近似解之表达式代入前述问题之弱形式,并取检验函数,可得

写成矩阵形式,便成为常见之有限元方程 由于历史原因,通常采用固体力学中的习惯命名:[公式] 为刚度矩阵,[公式] 为自由度向量,[公式] 为载荷向量。 由于空间已分片离散,上面的有限元方程只在各单元内部成立。为了求解方便,通常我们将所有单元的有限元方程连立起来求解,于是需要将各单元之刚度矩阵组装成总体刚度矩阵。 求解器MATLAB 代码及简释 本求解器之十行代码如下: function u=fem(nds,els,bcs) nnd=size(nds,1); u=zeros(nnd,1); K=zeros(nnd,nnd); f=zeros(nnd,1); for j=1:size(els,1) K(els(j,:),els(j,:))=K(els(j,:),els(j,:))+stima(nds(els(j,:),:)); f(els(j,:))=f(els(j,:))+ones(3,1)*det([1,1,1;nds(els(j,:),:)'])/ 6;

运用MATLAB对桁架单元进行有限元分析

np=3; ne=2; nload=1; nb=4; nu=0; % np为节点数,ne为单元数,nload为外力数,nb为约束数,nu 为泊松比 np2=np*2; np3=np2-nb; % np2为不受约束时自由度是节点的两倍 ,np3为不受约束的节点自由度个数 K=sym(zeros(np2,np2)); % 定义受整体刚度空数组 F=sym(zeros(np2,1)); % 定义节点外力空数组 KK=sym(zeros(np3,np3)); % 预置自由度刚度空数组 FF=sym(zeros(np3,1)); % 预置自由度外力空数组 syms h A P E L % 定义未知正常数为变量 xy=[0,h;sqrt(3)*h,0;0,0]; % 节点横纵坐标数组 AA=[A;A]; % 单元杆件面积 load=[2,2,-P]; % 荷载数组 bound=[1,1,0;1,2,0;3,1,0;3,2,0]; % 边界条件数组 IJ=[1,2;3,2]; % 各杆单元节点编号

DW=zeros(1,4); for ie =1:ne ip=IJ(ie,1); jp=IJ(ie,2); DW(1)=ip*2-1; DW(2)=ip*2; DW(3)=jp*2-1; DW(4)=jp*2; % 给单元节点横纵方向编号 x1=xy(ip,1); x2=xy(jp,1); y1=xy(ip,2); y2=xy(jp,2); % 杆单元节点x,y坐标 L=sqrt((x2-x1)^2+(y2-y1)^2); % 杆单元长度 c=(x2-x1)/L; s=(y2-y1)/L; % c为余弦, s为正弦 T=[c,s,0,0;-s,c,0,0;0,0,c,s;0,0,-s,c]; % 坐标转换矩阵(将局部坐标转换为整体坐标) A1=AA(ie); ke=[E*A1/L,0,-E*A1/L,0;0,0,0,0;-E*A1/L,0,E*A1/L,0; 0,0,0,0]; k=T.'*ke*T; % 将局部坐标下单元刚度转换为整体坐标下单元刚度 (转置后面加.可以去掉结果中的虚数) for i =1:4 i1=DW(i); for j =1:4 j1=DW(j); K(i1,j1)=K(i1,j1)+k(i,j); % 集成整体刚度矩阵

相关文档
相关文档 最新文档