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有限元分析与应用有限元分析(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有限元分析实例
1.物理现象:这个对工程师来说是直观的物理现象和物理量,温度多少度,载荷是多大等等。
通常来说,用户界面中呈现的、用户对工程问题进行设置时输入的都是此类信息。
2.数学方程:将物理现象翻译成相应的数学方程,例如流体对应的是NS方程,传热对应的是传热方程等等;大部分描述这些现象的方程在空间上都是偏微分方程,偶尔也有ODE(如粒子轨迹、化学反应等)。
在这个层面,软件把物理现象“翻译”为以解析式表示的数学模型。
3.数值模型:在定义了数学模型,并执行了网格剖分后,商业软件会将数学模型离散化,利用有限元方法、边界元法、有限差分法、不连续伽辽金法等方法生成数值模型。
软件会组装并计算方程组雅可比矩阵,并利用求解器求解方程组。
这个层面的计算通常是隐藏在后台的,用户只能通过一些求解器的参数来干预求解。
有限元是一种数值求解偏微分方程的方法。
基本过程大致是设置形函数,离散,形成求解矩阵,数值解矩阵,后处理之类的。
MATLAB要把这些过程均自己实现,不过在数值求解矩阵时可以调用已有函数。
可以理解为MATLAB是一个通用的计算器,当然它的功能远不止如此。
而ANSYS之类的叫做通用有限元软件,针对不同行业已经将上述过程封装,前后处理也比较漂亮,甚至不太了解有限元理论的人也能算些简单的东西,当然结果可靠性又另说了。
比较两者,ANSYS之类的用起来容易得多,但灵活性不如MATLAB。
MATLAB用起来很困难,也有人做了一些模块,但大多数只能解决一些相对简单的问题。
对于大多数工程问题,以及某些领域的物理问题,一般都用通用有限元软件,这些软件还能添加一些函数块,用以解决一些需要额外设置的东西。
但是对于非常特殊的问题,以及一般性方程的有限元解,那只能用MATLAB或C,Fortran之类的了。
桁架结构的有限元分析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 代码及简释本求解器之十行代码如下: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有限元法分析
一维梁的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 提供的 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有限元分析与应用》是2004年4月清华大学出版社出版的图书,作者是卡坦,译者是韩来彬。
内容简介:
《MATLAB有限元分析与应用》特别强调对MATLAB的交互应用,书中的每个示例都以交互的方式求解,使读者很容易就能把MATLAB用于有限分析和应用。
另外,《MATLAB有限元分析与应用》还提供了大量免费资源。
《MATLAB有限元分析与应用》采用当今在工程和工程教育方面非常流行的数学软件MATLAB来进行有限元的分析和应用。
《MATLAB有限元分析与应用》由简单到复杂,循序渐进地介绍了各种有限元及其分析与应用方法。
书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题,具有很高的工程应用价值。
matlab有限元分析实例
matlab有限元分析实例有限元是一种数值求解偏微分方程的方法。
基本过程大致是设置形函数,离散,形成求解矩阵,数值解矩阵,后处理之类的。
MATLAB要把这些过程均自己实现,不过在数值求解矩阵时可以调用已有函数。
可以理解为MATLAB是一个通用的计算器,当然它的功能远不止如此。
而ANSYS之类的叫做通用有限元软件,针对不同行业已经将上述过程封装,前后处理也比较漂亮,甚至不太了解有限元理论的人也能算些简单的东西,当然结果可靠性又另说了。
比较两者,ANSYS之类的用起来容易得多,但灵活性不如MATLAB。
MATLAB用起来很困难,也有人做了一些模块,但大多数只能解决一些相对简单的问题。
对于大多数工程问题,以及某些领域的物理问题,一般都用通用有限元软件,这些软件还能添加一些函数块,用以解决一些需要额外设置的东西。
但是对于非常特殊的问题,以及一般性方程的有限元解,那只能用MATLAB或C,Fortran之类的了。
matlab的优势是提供一个矩阵化的编程语言,matlab采用COO格式存储稀疏矩阵,在组装全局刚度阵的时候,可以由下标直接组装访问元素,处理单个单元的时候,对于二维模型,需要至少二重循环,加上对单元的遍历,至少循环层数较多。
众所周知,matlab在处理循环时速度较慢,特别是循环嵌套,速度更慢。
采用parfor,也不成熟。
优势就是在矩阵求解,雅可比矩阵与行列式求解的时候,有函数,又快又香。
对于C++,一般采用CSR格式存储稀疏阵,矩阵的行列,非零元,都要注意其存储位置。
编程求解比较费劲。
在处理单个单元的的时候行列式,小矩阵运算,都要自己编程,这个时候逻辑一定要清晰,代码相关的注释一定要写好。
否则编着编者,就把自己绕进去。
采用C++的福音是做并行时更容易控制,采用openmpi并行处理灵活。
在线性方程求解时,intel MKL是另外一个福音,intel的处理器还是intel理解,真的超级快,有机会你可以试一下,体会一下CPU 所有线程飙到100%的快感。
matlab有限元分析实例
为了使用matlab进行有限元分析,我们必须首先学习MATLAB 的基本操作以及使用MATLAB进行有限元分析的基本操作。
复习:最后一课分析弹簧系统,得出系统刚度矩阵matlab有限元分析20140226 matlab有限元分析20140226为了使用matlab进行有限元分析,我们首先必须学习MATLAB的基本操作,还要学习使用matlab进行有限元分析的基本操作。
1.回顾:最后一堂课分析了弹簧系统并推导出系统刚度矩阵。
2.为了使用matlab进行有限元分析,我们应该首先学习MATLAB的基本操作,然后再学习使用MATLAB进行有限元分析的基本操作。
Matlab有限元分析20140226为了使用Matlab进行有限元分析,首先必须学习MATLAB的基本操作以及使用MATLAB进行有限元分析的基本操作。
1.回顾:在上一课中对弹簧系统进行了分析,并得出了系统刚度矩阵。
2. MATLAB有限元分析的基本运算单位划分(选择哪个单位和划分多少个元素)分为两部分:单位选择和单位数回顾:最后一课分析了运算的基础MATLAB弹簧系统有限元分析matlab有限元分析20140226为了使用matlab进行有限元分析,我们首先应该学习MATLAB的基本操作,也要学会使用matlab进行有限元分析的基本操作。
1.回顾:最后一堂课分析了弹簧系统并推导出系统刚度矩阵。
2. Matlab导出系统刚度矩阵。
Matlab有限元分析的操作基础Matlab有限元分析20140226为了使用MATLAB进行有限元分析,首先必须学习MATLAB的基本操作,还必须学会使用Matlab 进行有限元分析的基本操作。
matlab 有限元 编程
MATLAB 是一种高级编程语言和交互式环境,常用于数值计算、数据分析、算法开发、数据可视化以及应用开发等。
在有限元分析(FEA)中,MATLAB 可以用来编写和运行有限元程序。
以下是一个简单的 MATLAB 有限元分析编程示例:这个例子假设你正在解决一个一维拉普拉斯方程,其形式为 -d^2u/dx^2 = f,在区间 [0, 1] 上。
```matlab% 参数定义L = 1; % 长度h = 0.01; % 步长N = L/h; % 元素数量x = linspace(0, L, N+1); % x 向量% 创建有限元网格nodes = x(1:N+1);elements = [nodes(2:end-1) nodes(1:end-1) nodes(2:end)]; % 定义载荷向量 ff = sin(pi*x);% 定义刚度矩阵 A 和载荷向量 FA = zeros(N, N);F = zeros(N, 1);for i = 1:Nfor j = 1:i-1A(i, j) = A(j, i) = h^2/(6*(x(i)-x(j))^2);endF(i) = -f(i)*h/2;end% 使用 MATLAB 的线性方程求解器u = A\F;% 绘制结果plot(x, u);xlabel('x');ylabel('u');title('有限元解');```这个代码首先定义了问题的参数,然后创建了一个有限元网格。
然后,它定义了刚度矩阵 A 和载荷向量 F。
最后,它使用 MATLAB 的线性方程求解器来求解方程,并绘制结果。
请注意,这只是一个非常简单的例子。
在实际的有限元分析中,问题可能会更复杂,并且需要更多的编程工作。
有限元的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”对话框。
Matlab有限元分析操作基础
Matlab有限元分析操作基础Matlab 有限元分析20140226为了用Matlab 进行有限元分析,首先要学会Matlab 基本操作,还要学会使用Matlab 进行有限元分析的基本操作。
1. 复习:上节课分析了弹簧系统x 推导了系统刚度矩阵1122121200k 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是多大矩阵?今天的系统刚度矩阵是什么?因为11221212k kk kk k k k----+所以10001000200200 100200300----c) K=SpringAssemble(K,k1,1,2) function y = SpringAssemble(K,k,i,j) 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); 1100100100100k -??=??-??10010001001000000K -??=-??K=SpringAssemble(K,k2,2,3) 1200200200200k -??10010001003002000200200K -??=---??1001000100010010030020002002000200200100200300----≠----步骤四:引入边界条件,消除冗余方程>>k=K(2:3,2:3)%构造不含冗余的方程>>f=[0;15]%构造外力列阵步骤五:解方程引例:已知1212u 31u u u +=??-=?,求 12u u 和解:类似求解KU=F ,输入下列Matlab 命令:>> K=[1 1;1,-1]>> F=[3;1]>> U=inv(K)*F>> U=K \F(继续弹簧系统求解)>>u=k \f %使用高斯消去法求解>>U=[0 ; u]%构造原方程组>>F=K*U %求出所有外力,含多余计算步骤六:后处理、扩展计算>>u1=[0;U(2)]%构造单元位移>>f1=SpringElementForces(k1,u1)%求单元1内力>>u2=[U(2) ; U(3)]%构造单元2位移>>f2=SpringElementForces(k2,u2)%求单元2内力4. 总结clck1=SpringElementStiffness(100)%创建单元刚度矩阵 1 k2=SpringElementStiffness(200)%创建单元刚度矩阵 2 K=zeros(3,3)%创建空白整体刚度矩阵K=SpringAssemble(K,k1,1,2)%按节点装入单元矩阵 1 K=SpringAssemble(K,k2,2,3)%按节点装入单元矩阵2 k=K(2:3,2:3)%构造不含冗余的方程f=[0;15]%构造外力列阵u=k\f%使用高斯消去法求解U=[0 ; u]%构造系统节点位移列阵F=K*U%求出所有外力,含多余计算u1=[0;U(2)]%构造单元位移f1=SpringElementForces(k1,u1)%求单元1内力u2=[U(2) ; U(3)]%构造单元2位移f2=SpringElementForces(k2,u2)%求单元2内力5. 练习1 Danyi 132 dan 34 3dan 35 4dan 35 dan5 54 dan6 42。
matlab有限元法圆柱绕流
matlab有限元法圆柱绕流分析,需要使用到MATLAB的有限元分析工具箱(FEATool)。
以下是一个简单的步骤:
1.建立模型:首先,使用FEATool创建圆柱绕流的模型。
这包括
定义流体区域、边界条件、初始条件等。
2.划分网格:使用FEATool对模型进行网格划分,以便进行有限
元分析。
3.设置材料属性:为流体和圆柱设置相应的物理属性,如密度、
粘度等。
4.设置边界条件和载荷:定义流体的入口和出口速度,以及任何
作用在圆柱上的力或力矩。
5.运行分析:使用FEATool进行有限元分析,这可能涉及到求解
流场的动量方程和连续性方程等。
6.后处理:使用FEA Tool的结果可视化功能,查看流场的速度、
压力等分布,以及圆柱受到的力或力矩等。
请注意,这只是一个非常基础的概述。
实际操作中,可能需要考虑更多的因素,如湍流模型的选择、边界条件的详细设置等。
另外,由于有限元分析是一个计算密集型任务,可能需要高性能计算资源。
matlab有限元分析一维弹簧系统
图 一维弹簧系统 提示:[]e K K K K K -?? =??-?? 函数脚本1 function y = SpringElementStiffness(k) y = [k, -k; -k, k]; end 函数脚本2 function y = SpringAssemble(K, k, i, j) 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; end 函数脚本3 function y = SpringElementForces(k,u) y = k*u; end 命令程序 >> k1=SpringElementStiffness(100); >> k2=SpringElementStiffness(200); >> n=3; >> K=zeros(n,n); >> K=SpringAssemble(K,k1,1,2); >> K=SpringAssemble(K,k2,2,3) >> U=zeros(2,1); >> K1 = K(2:3,2:3); >> U = K1\F; >> U=[0;U];
文档之家的所有文档均为用户上传分享文档之家仅负责分类整理如有任何问题可通过上方投诉通道反馈
matlab有 限 元 分 析 一 维 弹 簧 系 统
作业:编程求解(2019.4.28) 一维弹簧系统,已知K 1=100lb/in ,K 2=200lb/in ,F=300lb ,求: ① 整体刚度矩阵;②右端点的位移;③约束反力;④弹簧 的内力
有限元的matlab编程
本程序的网架位移求解函数附在主程序后面,主程序运行时调用该函数。
完整ppt
14
几何建模 定义荷载
用自定义输入
加下划线的为单元编号
完整ppt
5
形成等效荷载列阵
f=[0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a;0;0;0;a];%每个节点两个自由度,a 为之前输入的节点力
集成总刚:
获取单元两端节点坐标
xi = Node( Element( ie, 1 ), 1 ) ;%ie为单元号,以下相同 yi = Node( Element( ie, 1 ), 2 ) ; xj = Node( Element( ie, 2 ), 1 ) ; yj = Node( Element( ie, 2 ), 2 ) ;
节点号:',num2str(Element(ie,2)),' 轴力:',num2str(nodef(1))] ) ;
完整ppt
12
end
例二:网架
完整ppt
13
思路分析
几点说明
网架是由多根杆件按照一定的网格形式通过节点连结而成的空间结构。构成 网架的基本单元有三角锥,三棱体,正方体,截头四角锥等。鉴于网架的形 式较多,本程序提供一种通用的网架输入方法,但录入较为繁琐,同时提供 一种正放四角锥网架的简易输入方法作为典型。
有限元编程示例
完整ppt
1
例一:桁架
题目描述:
如下图所示的平面桁架,杆件长度、弹性模量、截 面积以及所受节点力P的大小可以自行定义。求节 点位移及杆件轴力。
《有限元基础教程》_【MATLAB算例】4.8.1(1) 基于4节点四面体单元的空间块体分析(Tetrahedron3D4Node)
【MATLAB 算例】4.8.1(1) 基于4节点四面体单元的空间块体分析(Tetrahedron3D4Node)如图4-22所示的一个块体,在右端面上端点受集中力F 作用。
基于MATLAB 平台,计算各个节点位移、支反力以及单元的应力。
取相关参数为:10110Pa,=0.25E μ=⨯,5=110N F ⨯。
图4-22 一个空间块体的分析解答:对该问题进行有限元分析的过程如下。
(1)结构的离散化与编号将结构离散为5个4节点四面体单元,单元编号及节点编号和坐标如图4-22所示,连接关系见表4-8,节点的坐标见表4-9。
表4-8 单元连接关系单元号 节点号 1 2 3 4 51 42 6 1 43 7 6 7 5 1 6 7 84 1 4 6 7表4-9 节点的坐标节点节点坐标/mxyz 1 2 3 4 5 6 7 80 0 0 0.2 0 0 0 0.8 0 0.2 0.8 0 0 0 0.6 0.2 0 0.6 0 0.8 0.6 0.20.80.6节点位移列阵[]111222888 Tu v w u v w u v w =q (4-190)节点外载列阵34780 0 0 0 TT T T T⎡⎤=⎣⎦F F F F F(4-191)其中34785000 00110N ⎡⎤⎡⎤⎢⎥⎢⎥====⎢⎥⎢⎥⎢⎥⎢⎥-⨯⎣⎦⎣⎦F F F F约束的支反力列阵12560000TTT T T ⎡⎤=⎣⎦R R R R R(4-192其中1256112255661256 x x x x y y y y z z z z R R R R R R R R R R R R ⎡⎤⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥====⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦⎣⎦R R R R总的节点载荷列阵12345678 TT T T T T T T T⎡⎤=+==⎣⎦P F R R R R F F R R F F (4-193)(2)计算各单元的刚度矩阵(以国际标准单位)首先在MA TLAB 环境下,输入弹性模量E 、泊松比NU ,然后针对单元1和单元2,分别5次调用函数Tetrahedron3D4Node_Stiffness ,就可以得到单元的刚度矩阵k1(6×6) ~ k5(6×6)。
悬臂梁MATLAB有限元算例注释
用有限元法对悬臂梁分析的算例算例:如下图所示的悬臂梁,受均布载荷q =1N /mm 2作用。
E =2.1×105N /mm 2, μ=0.3厚度h =10mm 。
现用有限元法分析其位移及应力。
梁可视为平面应力状态,先按图示尺寸划分为均匀的三角形网格,共有8×10=80个单元,5×ll =55个节点,坐标轴以及单元与节点的编号如图。
将均布载荷分配到各相应节点上,把有约束的节点5l 、52、53、54、55视作固定铰链,建立如图所示的离散化计算模型。
程序计算框图:(续左)(接右)开 始 输入材料参数 计算具有代表性的单元刚阵 K<=0 将各单元刚阵按整体编号集成到整体刚阵 处理根部约束,修改【K 】【Q 】 求解[K][δ]=[Q] 整理[δ] 并画图计算单元应力,并输出结束程序中的函数功能介绍及源代码1. LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym)――该函数用于计算平面应力情况下弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)、第二个节点坐标为(xj,yj)、第三个节点坐标为(xm,ym)时的线性三角形元的单元刚度矩阵.该函数返回6×6的单位刚度矩阵k.2. LinearTriangleAssemble(K,k,i,j,m)――该函数将连接节点i,j,m的线性三角形元的单元刚度矩阵k集成到整体刚度矩阵K。
每集成一个单元,该函数都将返回2N×2N的整体刚度矩阵K.3. LinearTriangleElementStresses(E,NU,t,xi,yi,xj,yj,xm,ym,u)-- 该函数计算在平面应力情况下弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)第二个节点坐标为(xj,yj)、第三个节点坐标为(xm,ym)以及单元位移矢量为u时的单元应力。
2012年-有限元作业-matlab编程实现有限元求解简单结构位移及应力
2012年有限元作业对于三角形单元有: 1、形函数⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎣⎡⎥⎦⎤=⎩⎨⎧⎭⎬⎫k k j j i i k jik j i k j i k j ik j ik j i v u v u v u c c c b b b a a a c c c b b b a a a y x y x A v u 00000000000000001000000121 其中:j m m i m m j j i y x y x y x y x a -==m j mji y y y y b -=-=11j m m j i x x x x c -==11 kkj ji iy x y x y x A 1112=相乘后得:()])()()[(21,k k k k j j j j i i i i u y c x b a u y c x b a u y c x b a A y x u ++++++++=()])()()[(21,k k k k j j j j i i i i v y c x b a v y c x b a v y c x b a Ay x v ++++++++=也可写成:()()kk j j i i k k j j i i v N v N v N y x v u N u N u N y x u ++=++=,,简写为:][}{q N v u =⎩⎨⎧⎭⎬⎫其中}}{{Tk kj j i iv u v u v uq =Ay c x b a N A y c x b a N A y c x b a N k k k i j j j i i i i i 2/)(2/)(2/)(++=++=++= 2、应变有弹性力学知道:xvy u y v x u xy y x ∂∂+∂∂=∂∂=∂∂==γεε,,,可得 }{}{}{q B v u v u v u c b c b c b c c c b b b A u b u b u b v c v c v c v c v c v c u b u b u b A v u x yy xy v x u y v x u k k j j i i k kjjiik j i k j ik k j j i i k k j j i i k k j j i i k k j j i i =⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤=⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤+++++++++=⎩⎨⎧⎭⎬⎫⎪⎪⎪⎩⎪⎪⎪⎨⎧⎪⎪⎪⎭⎪⎪⎪⎬⎫∂∂∂∂∂∂∂∂=⎪⎪⎪⎩⎪⎪⎪⎨⎧⎪⎪⎪⎪⎪⎪⎪⎭⎫∂∂+∂∂∂∂∂∂=000000212100ε3、应力}{][}{[][]{}q B D D ==εσ[]⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=2100010112μμμμE D 4、刚度矩阵[]()⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-+-+-+-+-=s r s r sr s r s r s r s r s r rs b b c c cb bc b c b b c c b b A Et K 21212121142μμμμμμμ 题目:如图所示是一平面梁。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 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有限元分析与应用》由简单到复杂,循序渐进地介绍了各种有限元及其分析与应用方法。
书中提供了大量取自机械工程、土木工程、航空航天工程和材料科学的示例和习题,具有很高的工程应用价值。