matlab有限元分析实例

合集下载

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

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

有限元大作业程序设计学校:天津大学院系:建筑工程与力学学院专业:01级工程力学姓名:刘秀学号:\\\\\\\\\\\指导老师:连续体平面问题的有限元程序分析[题目]:如图所示的正方形薄板四周受均匀载荷的作用,该结构在边界上受正向分布压力,m kNp 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程序的原理如框图:(1)主要变量:ID:问题类型码,ID=1时为平面应力问题,ID=2时为平面应变问题N_NODE:节点个数N_LOAD:节点载荷个数N_DOF:自由度,N_DOF=N_NODE*2(平面问题)N_ELE:单元个数N_BAND:矩阵半带宽N_BC:有约束的节点个数PE:弹性模量PR:泊松比PT:厚度LJK_ELE(I,3):单元节点编号数组,LJK_ELE(I,1),LJK_ELE(I,2),LJK_ELE(I,3)分别放单元I的三个节点的整体编号X(N_NODE), Y(N_NODE):节点坐标数组,X(I),Y(I)分别存放节点I的x,y 坐标值P_LJK(N_BC,3):节点载荷数组,P_LJK(I,1)表示第I个作用有节点载荷的节点的编号,P_LJK(I,2),P_LJK(I,3)分别为该节点沿x,y方向的节点载荷数值AK(N_DOF,N_BAND):整体刚度矩阵AKE(6,6):单元刚度矩阵BB(3,6):位移……应变转换矩阵(三节点单元的几何矩阵)DD(3,3):弹性矩阵SS(3,6);应力矩阵RESULT_N(N_NOF):节点载荷数组,存放节点载荷向量,解方程后该矩阵存放节点位移DISP_E(6)::单元的节点位移向量STS_ELE(N_ELE,3):单元的应力分量STS_ND(N_NODE,3):节点的应力分量(2)子程序说明:READ_IN:读入数据BAND_K:形成半带宽的整体刚度矩阵FORM_KE:计算单元刚度矩阵FORM_P:计算节点载荷CAL_AREA:计算单元面积DO_BC:处理边界条件CLA_DD:计算单元弹性矩阵SOLVE:计算节点位移CLA_BB:计算单元位移……应变关系矩阵CAL_STS:计算单元和节点应力(3)文件管理:源程序文件:chengxu.for程序需读入的数据文件:BASIC.IN,NODE.IN,ELEMENT.IN(需要手工生成)程序输出的数据文件:DATA.OUT(4)数据文件格式:需读入的模型基本信息文件BASIC.IN的格式如下表需读入的节点信息文件NODE.IN的格式如下表需读入的单元信息文件ELEMENT.IN的格式如下表输出结果文件DATA.OUT格式如下表[算例原始数据和程序分析]:(1)模型基本信息文件BASIC.IN的数据为1,4,6,5,31.,0.,1.1,1,0,2,1,0,4,1,1,5,0,1,6,0,11,-0.5,-1.5,3.,-1.,-1,6,-0.5,-0.5(2)手工准备的节点信息文件NODE.IN的数据为1 0.0 2.02 0.0 1.03 1.0 1.04 0. 0.5 1.0 0.6 2.0 0.(3)手工准备的单元信息文件ELEMENT.IN的数据为1 2 3 3 0 0 0 0 1 1 1 1 0 12 4 5 5 0 0 0 0 1 1 1 1 0 25 3 2 2 0 0 0 0 1 1 1 1 0 33 5 6 6 0 0 0 0 1 1 1 1 04 (4)源程序文件chengxu.for为:PROGRAM FEM2DDIMENSION IJK_ELE(500,3),X(500),Y(500),IJK_U(50,3),P_IJK(50,3),&RESULT_N(500),AK(500,100)D IMENSION STS_ELE(500,3),STS_ND(500,3)OPEN(4,FILE='BASIC.IN')OPEN(5,FILE='NODE.IN')OPEN(6,FILE='ELEMENT.IN')OPEN(8,FILE='DATA.OUT')OPEN(9,FILE='FOR_POST.DAT')READ(4,*)ID,N_ELE,N_NODE,N_BC,N_LOADIF(ID.EQ.1)WRITE(8,20)IF(ID.EQ.2)WRITE(8,25)20 FORMAT(/5X,'=========PLANE STRESS PROBLEM========')25 FORMAT(/5X,'=========PLANE STRAIN PROBLEM========')CALL READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR,PT, & IJK_ELE,X,Y,IJK_U,P_IJK)CALL BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,& IJK_ELE,X,Y,PE,PR,PT,AK)CALL FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)CALL DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N)CALL SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N)CALL CAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, & STS_ELE,STS_ND)c to putout a data fileWRITE(9,70)REAL(N_NODE),REAL(N_ELE)70 FORMAT(2f9.4)WRITE(9,71)(X(I),Y(I),RESULT_N(2*I-1),RESULT_N(2*I),& STS_ND(I,1),STS_ND(I,2),STS_ND(I,3),I=1,N_NODE)71 FORMA T(7F9.4)WRITE(9,72)(REAL(IJK_ELE(I,1)),REAL(IJK_ELE(I,2)),&REAL(IJK_ELE(I,3)),REAL(IJK_ELE(I,3)),&STS_ELE(I,1),STS_ELE(I,2),STS_ELE(I,3),I=1, N_ELE)72 FORMAT(7f9.4)cCLOSE(4)CLOSE(5)CLOSE(6)CLOSE(8)CLOSE(9)E NDcc to get the original data in order to model the problemSUBROUTINE READ_IN(ID,N_ELE,N_NODE,N_BC,N_BAND,N_LOAD,PE,PR, &PT,IJK_ELE,X,Y,IJK_U,P_IJK)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),IJK_U(N_BC,3), & P_IJK(N_LOAD,3),NE_ANSYS(N_ELE,14)REAL ND_ANSYS(N_NODE,3)READ(4,*)PE,PR,PTREAD(4,*)((IJK_U(I,J),J=1,3),I=1,N_BC)READ(4,*)((P_IJK(I,J),J=1,3),I=1,N_LOAD)READ(5,*)((ND_ANSYS(I,J),J=1,3),I=1,N_NODE)READ(6,*)((NE_ANSYS(I,J),J=1,14),I=1,N_ELE)DO 10 I=1,N_NODEX(I)=ND_ANSYS(I,2)Y(I)=ND_ANSYS(I,3)10 CONTINUEDO 11 I=1,N_ELEDO 11 J=1,3IJK_ELE(I,J)=NE_ANSYS(I,J)11 CONTINUEN_BAND=0DO 20 IE=1,N_ELEDO 20 I=1,3DO 20 J=1,3IW=IABS(IJK_ELE(IE,I)-IJK_ELE(IE,J))IF(N_BAND.LT.IW)N_BAND=IW20 CONTINUEN_BAND=(N_BAND+1)*2IF(ID.EQ.1) THENELSEPE=PE/(1.0-PR*PR)PR=PR/(1.0-PR)END IFR ETURNENDcC to form the stiffness matrix of elementSUBROUTINE FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6),DD(3,3), & AKE(6,6), SS(6,6)CALL CAL_DD(PE,PR,DD)CALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 10 I=1,3DO 10 J=1,6SS(I,J)=0.0DO 10 K=1,310 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 20 I=1,6DO 20 J=1,6AKE(I,J)=0.0DO 20 K=1,320 AKE(I,J)=AKE(I,J)+SS(K,I)*BB(K,J)*AE*PTRETURNENDcc to form banded global stiffness matrixSUBROUTINE BAND_K(N_DOF,N_BAND,N_ELE,IE,N_NODE,IJK_ELE,X,Y,PE, & PR,PT,AK)DIMENSIONIJK_ELE(500,3),X(N_NODE),Y(N_NODE),AKE(6,6),AK(500,100)N_DOF=2*N_NODEDO 40 I=1,N_DOFDO 40 J=1,N_BAND40 AK(I,J)=0DO 50 IE=1,N_ELECALL FORM_KE(IE,N_NODE,N_ELE,IJK_ELE,X,Y,PE,PR,PT,AKE)DO 50 I=1,3DO 50 II=1,2IH=2*(I-1)+IIIDH=2*(IJK_ELE(IE,I)-1)+IIDO 50 J=1,3DO 50 JJ=1,2IL=2*(J-1)+JJIZL=2*(IJK_ELE(IE,J)-1)+JJIDL=IZL-IDH+1IF(IDL.LE.0) THENELSEAK(IDH,IDL)=AK(IDH,IDL)+AKE(IH,IL)END IF50 CONTINUERETURNENDcc to calculate the area of elementSUBROUTINE CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)XIJ=X(J)-X(I)YIJ=Y(J)-Y(I)XIK=X(K)-X(I)YIK=Y(K)-Y(I)AE=(XIJ*YIK-XIK*YIJ)/2.0RETURNENDcc to calculate the elastic matrix of elementSUBROUTINE CAL_DD(PE,PR,DD)DIMENSION DD(3,3)DO 10 I=1,3DO 10 J=1,310 DD(I,J)=0.0DD(1,1)=PE/(1.0-PR*PR)DD(1,2)=PE*PR/(1.0-PR*PR)DD(2,1)=DD(1,2)DD(2,2)=DD(1,1)DD(3,3)=PE/((1.0+PR)*2.0)RETURNENDcc to calculate the strain-displacement matrix of elementSUBROUTINE CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB) DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),BB(3,6)I=IJK_ELE(IE,1)J=IJK_ELE(IE,2)K=IJK_ELE(IE,3)DO 10 II=1,3DO 10 JJ=1,310 BB(II,JJ)=0.0BB(1,1)=Y(J)-Y(K)BB(1,3)=Y(K)-Y(I)BB(1,5)=Y(I)-Y(J)BB(2,2)=X(K)-X(J)BB(2,4)=X(I)-X(K)BB(2,6)=X(J)-X(I)BB(3,1)=BB(2,2)BB(3,2)=BB(1,1)BB(3,3)=BB(2,4)BB(3,4)=BB(1,3)BB(3,5)=BB(2,6)BB(3,6)=BB(1,5)CALL CAL_AREA(IE,N_NODE,IJK_ELE,X,Y,AE)DO 20 I1=1,3DO 20 J1=1,620 BB(I1,J1)=BB(I1,J1)/(2.0*AE)RETURNENDcc to form the global load matrixSUBROUTINE FORM_P(N_ELE,N_NODE,N_LOAD,N_DOF,IJK_ELE,X,Y,P_IJK, & RESULT_N)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),P_IJK(N_LOAD,3), & RESULT_N(N_DOF)DO 10 I=1,N_DOF10 RESULT_N(I)=0.0DO 20 I=1,N_LOADII=P_IJK(I,1)RESULT_N(2*II-1)=P_IJK(I,2)20 RESULT_N(2*II)=P_IJK(I,3)RETURNENDcc to deal with BC(u) (here only for fixed displacement) using "1-0" method SUBROUTINE DO_BC(N_BC,N_BAND,N_DOF,IJK_U,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),IJK_U(N_BC,3),AK(500,100)DO 30 I=1,N_BCIR=IJK_U(I,1)DO 30 J=2,3IF(IJK_U(I,J).EQ.0)THENELSEII=2*IR+J-3AK(II,1)=1.0RESULT_N(II)=0.0DO 10 JJ=2,N_BAND10 AK(II,JJ)=0.0DO 20 JJ=2,II20 AK(II-JJ+1,JJ)=0.0END IF30 CONTINUERETURNENDcc to solve the banded FEM equation by GAUSS eliminationSUBROUTINE SOLVE(N_NODE,N_DOF,N_BAND,AK,RESULT_N) DIMENSION RESULT_N(N_DOF),AK(500,100)DO 20 K=1,N_DOF-1IF(N_DOF.GT.K+N_BAND-1)IM=K+N_BAND-1IF(N_DOF.LE.K+N_BAND-1)IM=N_DOFDO 20 I=K+1,IML=I-K+1C=AK(K,L)/AK(K,1)IW=N_BAND-L+1DO 10 J=1,IWM=J+I-K10 AK(I,J)=AK(I,J)-C*AK(K,M)20 RESULT_N(I)=RESULT_N(I)-C*RESULT_N(K)RESULT_N(N_DOF)=RESULT_N(N_DOF)/AK(N_DOF,1)DO 40 I1=1,N_DOF-1I=N_DOF-I1IF(N_BAND.GT.N_DOF-I-1)JQ=N_DOF-I+1IF(N_BAND.LE.N_DOF-I-1)JQ=N_BANDDO 30 J=2,JQK=J+I-130 RESULT_N(I)=RESULT_N(I)-AK(I,J)*RESULT_N(K)40 RESULT_N(I)=RESULT_N(I)/AK(I,1)WRITE(8,50)50 FORMAT(/12X,'* * * * * RESULTS BY FEM2D * * * * *',//8X,&'--DISPLACEMENT OF NODE--'//5X,'NODE NO',8X,'X-DISP',8X,'Y-DISP') DO 60 I=1,N_NODE60 WRITE(8,70) I,RESULT_N(2*I-1),RESULT_N(2*I)70 FORMAT(8X,I5,7X,2E15.6)RETURNENDcc calculate the stress components of element and nodeSUBROUTINECAL_STS(N_ELE,N_NODE,N_DOF,PE,PR,IJK_ELE,X,Y,RESULT_N, &STS_ELE,STS_ND)DIMENSION IJK_ELE(500,3),X(N_NODE),Y(N_NODE),DD(3,3),BB(3,6), &SS(3,6),RESULT_N(N_DOF),DISP_E(6)DIMENSION STS_ELE(500,3),STS_ND(500,3)WRITE(8,10)10 FORMAT(//8X,'--STRESSES OF ELEMENT--')CALL CAL_DD(PE,PR,DD)DO 50 IE=1,N_ELECALL CAL_BB(IE,N_NODE,N_ELE,IJK_ELE,X,Y,AE,BB)DO 20 I=1,3DO 20 J=1,6SS(I,J)=0.0DO 20 K=1,320 SS(I,J)=SS(I,J)+DD(I,K)*BB(K,J)DO 30 I=1,3DO 30 J=1,2IH=2*(I-1)+JIW=2*(IJK_ELE(IE,I)-1)+J30 DISP_E(IH)=RESULT_N(IW)STX=0STY=0TXY=0DO 40 J=1,6STX=STX+SS(1,J)*DISP_E(J)STY=STY+SS(2,J)*DISP_E(J)40 TXY=TXY+SS(3,J)*DISP_E(J)STS_ELE(IE,1)=STXSTS_ELE(IE,2)=STYSTS_ELE(IE,3)=TXY50 WRITE(8,60)IE,STX,STY,TXY60 FORMAT(1X,'ELEMENT NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)c the following part is to calculate stress components of nodeWRITE(8,55)55 FORMAT(//8X,'--STRESSES OF NODE--')DO 90 I=1,N_NODEA=0.B=0.C=0.II=0DO 70 K=1,N_ELEDO 70 J=1,3IF(IJK_ELE(K,J).EQ.I) THENII=II+1A=A+STS_ELE(K,1)B=B+STS_ELE(K,2)C=C+STS_ELE(K,3)END IF70 CONTINUESTS_ND(I,1)=A/IISTS_ND(I,2)=B/IISTS_ND(I,3)=C/IIWRITE(8,75)I,STS_ND(I,1),STS_ND(I,2),STS_ND(I,3)75 FORMAT(1X,'NODE NO.=',I5/18X,'STX=',E12.6,5X,'STY=',&E12.6,2X,'TXY=',E12.6)90 CONTINUERETURNENDc FEM2D programm end[算例结果]:chengxu.for所输出的数据文件DATA.OUT数据内容如下:=========PLANE STRESS PROBLEM========* * * * * RESULTS BY FEM2D * * * * *--DISPLACEMENT OF NODE--NODE NO X-DISP Y-DISP1 .000000E+00 -.525275E+012 .000000E+00 -.225275E+013 -.108791E+01 -.137363E+014 .000000E+00 .000000E+005 -.824176E+00 .000000E+006 -.182418E+01 .000000E+00--STRESSES OF ELEMENT--ELEMENT NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00ELEMENT NO.= 2STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00ELEMENT NO.= 3STX=-.108791E+01 STY=-.137363E+01 TXY= .307692E+00ELEMENT NO.= 4STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00--STRESSES OF NODE--NODE NO.= 1STX=-.108791E+01 STY=-.300000E+01 TXY= .439560E+00NODE NO.= 2STX=-.100000E+01 STY=-.220879E+01 TXY= .249084E+00NODE NO.= 3STX=-.105861E+01 STY=-.191575E+01 TXY= .205128E+00NODE NO.= 4STX=-.824176E+00 STY=-.225275E+01 TXY= .000000E+00NODE NO.= 5STX=-.970696E+00 STY=-.166667E+01 TXY= .586081E-01NODE NO.= 6STX=-.100000E+01 STY=-.137363E+01 TXY=-.131868E+00[结论与体会]:通过本次的课程设计,我对有限元的概念有了更加深刻的理解,同时也弥补了平时学习是疏忽的地方,充实了有限元知识。

第三章MATLAB有限元分析与应用

第三章MATLAB有限元分析与应用

第三章MATLAB有限元分析与应用有限元分析(Finite Element Analysis, FEA)是一种工程计算方法,用于解决结构力学和流体力学等问题。

它将一个复杂的结构分割成多个简单的离散单元,通过建立数学模型和求解方程组,得到结构的力学、热力学和流体力学等性能参数。

MATLAB是一种功能强大的数学计算软件,具有直观的用户界面和丰富的工具箱,可以方便地进行有限元分析。

本章将介绍在MATLAB中进行有限元分析的基本步骤和方法,以及一些常见的应用例子。

首先,进行有限元分析需要将结构进行离散化。

常用的离散化方法有节点法和单元法。

节点法是将结构的几何形状划分为小的节点,并在节点上进行计算。

单元法是将结构划分为多个小的单元,并在每个单元内进行计算。

在MATLAB中,可以通过创建节点和单元的矩阵来描述结构和单元的关系。

例如,创建一个2D结构形式的节点矩阵:nodes = [0 0; 1 0; 0 1; 1 1];然后,通过创建描述节点连接关系的矩阵,来定义结构的单元:elements = [1 2 3; 2 4 3];这里的每一行代表一个单元,数字表示节点的编号。

接下来,需要定义材料的力学参数和边界条件。

材料的力学参数包括弹性模量、泊松比等。

边界条件包括支座约束和加载条件。

在MATLAB中,可以通过定义力学参数和边界条件的向量来描述。

例如,定义弹性模量和泊松比的向量:E=[200e9200e9];%弹性模量nu = [0.3 0.3]; % 泊松比定义支座约束的向量(1表示固定,0表示自由):constraints = [1 1; 0 0; 0 1; 0 1];定义加载条件的向量(包括点力和面力):最后,通过求解方程组得到结构的应力和位移等结果。

在MATLAB中,可以利用有限元分析工具箱中的函数进行计算。

例如,可以使用“assem”函数将节点和单元的信息组装成方程组,并使用“solveq”函数求解方程组。

matlab有限元分析实例

matlab有限元分析实例

1.物理现象:这个对工程师来说是直观的物理现象和物理量,温度多少度,载荷是多大等等。

通常来说,用户界面中呈现的、用户对工程问题进行设置时输入的都是此类信息。

2.数学方程:将物理现象翻译成相应的数学方程,例如流体对应的是NS方程,传热对应的是传热方程等等;大部分描述这些现象的方程在空间上都是偏微分方程,偶尔也有ODE(如粒子轨迹、化学反应等)。

在这个层面,软件把物理现象“翻译”为以解析式表示的数学模型。

3.数值模型:在定义了数学模型,并执行了网格剖分后,商业软件会将数学模型离散化,利用有限元方法、边界元法、有限差分法、不连续伽辽金法等方法生成数值模型。

软件会组装并计算方程组雅可比矩阵,并利用求解器求解方程组。

这个层面的计算通常是隐藏在后台的,用户只能通过一些求解器的参数来干预求解。

有限元是一种数值求解偏微分方程的方法。

基本过程大致是设置形函数,离散,形成求解矩阵,数值解矩阵,后处理之类的。

MATLAB要把这些过程均自己实现,不过在数值求解矩阵时可以调用已有函数。

可以理解为MATLAB是一个通用的计算器,当然它的功能远不止如此。

而ANSYS之类的叫做通用有限元软件,针对不同行业已经将上述过程封装,前后处理也比较漂亮,甚至不太了解有限元理论的人也能算些简单的东西,当然结果可靠性又另说了。

比较两者,ANSYS之类的用起来容易得多,但灵活性不如MATLAB。

MATLAB用起来很困难,也有人做了一些模块,但大多数只能解决一些相对简单的问题。

对于大多数工程问题,以及某些领域的物理问题,一般都用通用有限元软件,这些软件还能添加一些函数块,用以解决一些需要额外设置的东西。

但是对于非常特殊的问题,以及一般性方程的有限元解,那只能用MATLAB或C,Fortran之类的了。

桁架结构的有限元分析MATLAB

桁架结构的有限元分析MATLAB

桁架结构的有限元分析MATLAB桁架结构是一种由直杆或斜杆连接而成的稳定结构,在工程应用中较为常见。

有限元分析(Finite Element Analysis,FEA)是一种利用数值方法解决结构力学问题的工具。

本文将介绍如何使用MATLAB进行桁架结构的有限元分析,并对其进行1200字以上的详细描述。

在进行桁架结构有限元分析前,需要先进行结构建模以及材料属性和加载条件的定义。

这些定义可以通过MATLAB命令行或者编写MATLAB脚本文件实现。

首先,我们需要定义桁架结构的节点和单元。

节点用于表示桁架结构的连接点,单元用于表示相邻节点之间的连接关系。

可以使用MATLAB中的矩阵表示节点和单元,如下所示:nodes = [x1, y1; x2, y2; ...; xn, yn];elements = [n1, n2; n3, n4; ...; nm, np];```其中,`nodes`是一个n行2列的矩阵,表示n个节点的坐标;`elements`是一个m行2列的矩阵,表示m个单元的节点连接关系。

接下来,我们需要定义材料属性和加载条件。

材料属性包括杨氏模量和截面面积等参数,加载条件包括节点的约束和外部加载。

可以使用MATLAB中的矩阵或者结构体表示材料属性和加载条件,如下所示:materials = struct('E', E1, 'A', A1; 'E', E2, 'A', A2; ...);constraints = [n1, d1x, d1y; ...; nm, dmx, dmy];loads = [n1, F1x, F1y; ...; nl, Flx, Fly];```其中,`materials`是一个结构体数组,每个结构体包含材料的杨氏模量(E)和截面面积(A);`constraints`是一个m行3列的矩阵,表示m个节点的约束,其中d1x和d1y分别表示节点的x方向和y方向位移约束;`loads`是一个l行3列的矩阵,表示l个节点的外部加载,其中F1x和F1y分别表示节点的x方向和y方向外部力。

matlab有限元分析实例

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;endfreends=setdiff(1:nnd,bcs);u(freends)=K(freends,freends)\f(freends);function stima=stima(vertices)ndim=size(vertices,2);J=[ones(1,ndim+1);vertices'];B=J\[zeros(1,ndim);eye(ndim)];stima=det(J).*B*B'/prod(1:ndim);输入数据格式:唯一的一个函数stima 用来计算各单元刚度矩阵。

一维梁的MATLAB有限元法分析

一维梁的MATLAB有限元法分析

一维梁的M A T L A B有限元法分析问题如下:梁A B在A和B两端固定,中间点表示为C,在中间区域承受均匀分布的载荷q,如图所示。

梁A B的抗弯刚度为E I。

1,使用R i t z法确定点C处的位移和弯矩,并讨论随着包含更多基本函数的准确性。

提示:根据偏转曲线的形状,可以选择基函数的形式,其中系数,并且应该由点A或B处的边界条件确定。

2,采用一维有限元法解决问题,并讨论网格越细时的准确性。

提示:使用1-D梁单元。

1R i t z法一维欧拉-伯努利梁的势能如下:设选择基函数,容易看出基函数满足边界条件设,代入势能表达式得到由于三角函数是正交函数系,所以得到令q=10N/m m,E=200000M P a,I=10000,L=200m m在M A T L A B中计算A k前十项得到A=-0.008400754770396-0.000320811945459-0.000049922673823 -0.000020050746591-0.000009258470170-0.000003960641302 -0.000001943426801-0.000001253171662-0.000000837688763 -0.000000513299113计算C点位移,使用1-10个试函数结果如下:0.0084007547703960.0084007547703960.0085006001180410.0085006001180410.0085191170583800.0085191170583800.0085230039119830.0085230039119830.0085246792895100.008524679289510计算C点弯矩,,使用1-10个试函数结果如下:0.8291212625436840.7024697829907620.746814316705690 0.7151514468174600.7379958063005830.723923419683592 0.7333220380018130.7254063205297550.732103122461494 0.727037063279377可以看到,位移收敛是很快的,弯矩收敛速度慢于位移。

matlab有限元分析实例

matlab有限元分析实例
唯一的一个函数stima 用来计算各单元刚度矩阵。为考虑进一步扩展之可能性,我们采用了等参单元的表达形式,引入了雅可比矩阵 J。事实上,该函数可直接用于三维空间四面体单元刚度矩阵之计算。
求解器之使用
可以使用文本文件存储输入数据,用MATLAB 提供的 load 命令载入。仅需下面三行代码即可完成读取输入数据,调用求解器求解,绘制结果云图诸个步骤。
ndim=size(vertices,2); J=[ones(1,ndim+1);vertices'];
B=J\[zeros(1,ndim);eye(ndim)]; stima=det(J).*B*B'/prod(1:ndim);
输入数据格式:
nds:节点列表矩阵。对于 [公式] 个节点,为一 [公式] 矩阵。第一列为 [公式] 坐标值,第二列为 [公式] 坐标值。
其中[公式] 为检验函数,[公式] 为一适当索伯列夫空间(在本例中也是一希尔伯特空间)。[公式] 为一双线性型,[公式] 为一线性型。其具体表达式为
[公式]
有限元空间离散
我们采用最简单的二维[公式] 单元离散 [公式]。[公式] 单元即三节点线性三角形单元,其插值基函数(即形函数) [公式] 为一次多项式。有限元之核心思想为:使用离散的函数空间 [公式] 来分片逼近连续的函数空间 [公式]。于是所求之近似解 [公式] 可写成基函数之线性组合
[公式]
其中[公式] 为待求系数,常称为自由度。将此近似解之表达式代入前述问题之弱形式,并取检验函数 [公式],可得
[公式]
写成矩阵形式,便成为常见之有限元方程
[公式]
由于历史原因,通常采用固体力学中的习惯命名:[公式] 为刚度矩阵,[公式] 为自由度向量,[公式] 为载荷向量。

matlab有限元分析实例

matlab有限元分析实例

有限元:
在数学中,有限元法是一种为求解偏微分方程边值问题近似解的数值技术。

求解时对整个问题区域进行分解,每个子区域都成为简单的部分,这种简单部分就称作有限元。

它通过变分方法,使得误差函数达到最小值并产生稳定解。

类比于连接多段微小直线逼近圆的思想,有限元法包含了一切可能的方法,这些方法将许多被称为有限元的小区域上的简单方程联系起来,并用其去估计更大区域上的复杂方程。

它将求解域看成是由许多称为有限元的小的互连子域组成,对每一单元假定一个合适的(较简单的)近似解,然后推导求解这个域总的满足条件(如结构的平衡条件),从而得到问题的解。

这个解不是准确解,而是近似解,因为实际问题被较简单的问题所代替。

MATLAB有限元分析与应用:
《MATLAB有限元分析与应用》是2004年4月清华大学出版社出版的图书,作者是卡坦,译者是韩来彬。

内容简介:
《MATLAB有限元分析与应用》特别强调对MATLAB的交互应用,书中的每个示例都以交互的方式求解,使读者很容易就能把MATLAB用于有限分析和应用。

另外,《MATLAB有限元分析与应用》还提供了大量免费资源。

《MATLAB有限元分析与应用》采用当今在工程和工程教育方面非常流行的数学软件MATLAB来进行有限元的分析和应用。

《MATLAB有限元分析与应用》由简单到复杂,循序渐进地介绍了各种有限元及其分析与应用方法。

书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题,具有很高的工程应用价值。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

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有限元分析与应用》由简单到复杂,循序渐进地介绍了各种有限元及其分析与应用方法。

书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题,具有很高的工程应用价值。

相关文档
最新文档