您的位置：360文档中心› 有限元MATLABWord版

有限元MATLABWord版

合集下载

有限元matlab仿真

有限元matlab仿真

有限元matlab仿真最近写有有限元作业，用matlab实现了整个的过程，顺便熟悉一下有限元的整体解题思路，感觉这确实是一种很好方法。

下面对书上的一道题做仿真吧，最后结果如下。

题目：如图所示薄板忽略自重，在力F的作用下，求各节点位移和各单元应力，已知假定该问题为平面应力问题。

这是一个正方形的薄板，上面两个顶点处被吊起来，在中心处受到一个向下的力。

现分成四个单元（下面的两个顶点编号为1,2，上面个两个为3,4，中间的5）执行程序得到如下结果1节点的位移:（u，v） = （-1.42857e-006, -7.14286e-006）2节点的位移:（u，v） = （1.42857e-006, -7.14286e-006）3节点的位移:（u，v） = （0, 0）4节点的位移:（u，v） = （0, 0）5节点的位移:（u，v） = （0, -8.57143e-006）第1单元:应力为（750000，3.75e+006，-2.25e+006）应变为（3.57143e-006，1.78571e-005，-2.14286e-005）第2单元:应力为（1.5e+006，-1.5e+006，6.98492e-010）应变为（7.14286e-006，-7.14286e-006，1.01644e-020）第3单元:应力为（750000，3.75e+006，2.25e+006）应变为（3.57143e-006，1.78571e-005，2.14286e-005）第4单元:应力为（0，9e+006，0）应变为（0，4.28571e-005，0）还可以实现各顶点移动的动画图，本来想做个gif，但发现getframe截不了屏，郁闷只能在这贴几张图了这个是不是就是ANSYS的有限元分析啊！代码部分如下，显示动画的代码麻烦点，就不贴了，从代码也可以看到，好的数据结构的重要性：main.m文件：clcclear allclose all% 数据choise = 1;switch choisecase 1% 作业第3题verts = [0 0;2 0;2 2;0 2;1 1;]*0.2;eleindex = [1 5 4;2 5 1;3 5 2;4 5 3];vertdisind = [1 1 1 1 0 0 0 0 0 1]'; R = [0 0 0 0 0 0 0 0 0 -30000]'; case 2% 书上例题verts = [0 2;0 1;1 1;0 0;1 0;2 0];eleindex = [5 2 4;6 3 5;2 5 3;3 1 2];endfigurepatch('Faces', eleindex, 'Vertices', verts, 'FaceColor', 'r'); % 绘制数据axis equalaxis off% 构造FEMhFEM.vertcount = size(verts, 1);hFEM.verts = verts;hFEM.elecount = size(eleindex, 1);hFEM.eleindex = eleindex;hFEM.mu = 0;hFEM.E = 210e9;hFEM.t = 0.01;hFEM = FEM(hFEM);% 求解有限元方程得各节点位移hFEM = FEMSolve(hFEM, vertdisind, R);for k = 1:hFEM.vertcountfprintf('%d节点的位移:（u，v）= （%g, %g）\n', k, hFEM.vertdis([2*k-1, 2*k]));end% 求各单元应力、应变for k = 1:hFEM.elecountind(1:2:5) = 2*hFEM.eleindex(k, :)-1;ind(2:2:6) = 2*hFEM.eleindex(k, :);stress = hFEM.eletruct(k).S * hFEM.vertdis(ind);strain = hFEM.eletruct(k).B * hFEM.vertdis(ind);fprintf('第%d单元:\n应力为（%g，%g，%g）\n', k, stress);fprintf('应变为（%g，%g，%g）\n', strain);endFEM.m文件：% 计算FEM分析的各种矩阵function hFEM = FEM(hFEM)% 处理每个单元for i = 1:hFEM.elecounthFEM.eletruct(i) = FEMele(hFEM, i);end% 计算合成刚度矩阵hFEM.K = FEMCalcCombiningK(hFEM);endfunction hFEMele = FEMele(hFEM, ind) hFEMele.ele = hFEM.verts(hFEM.eleindex(ind, :), :); hFEMele.B = FEMCalcB(hFEMele);hFEMele.D = FEMCalcD(hFEM);hFEMele.S = FEMCalcS(hFEMele);hFEMele.K = FEMCalcK(hFEM, hFEMele);end% 计算几何矩阵function B = FEMCalcB(hFEMele)ele = hFEMele.ele;A(:, 2:3) = ele;A(:, 1) = 1;invA = inv(A);B(1, 1:2:5) = invA(2, :);B(2, 2:2:6) = invA(3, :);B(3, 1:2:5) = invA(3, :);B(3, 2:2:6) = invA(2, :);end% 计算应力矩阵function D = FEMCalcD(hFEM)mu = hFEM.mu;E = hFEM.E;D = [1 mu 0;mu 1 0;0 0 (1-mu)/2]*E/(1-mu^2);end% 计算应力变换矩阵function S = FEMCalcS(hFEMele)S = hFEMele.D*hFEMele.B;end% 计算单元刚度矩阵function K = FEMCalcK(hFEM, hFEMele)triarea = CalcTriangleArea(hFEMele.ele);K = hFEM.t*triarea*hFEMele.B'*hFEMele.D*hFEMele.B; end% 计算三角形面积function area = CalcTriangleArea(ele)A(:, 2:3) = ele;A(:, 1) = 1;area = 0.5*det(A);end% 计算合成矩阵function K = FEMCalcCombiningK(hFEM) vertcount = hFEM.vertcount;elecount = hFEM.elecount;eleindex = hFEM.eleindex;eleK = zeros(6, 6, elecount);for k = 1:elecounteleK(:, :, k) = hFEM.eletruct(k).K;endK = zeros(vertcount*2);for i = 1:vertcountfor j = 1:vertcountfor k = 1:elecountindi = find(eleindex(k, :)==i);indj = find(eleindex(k, :)==j);if ~isempty(indi) && ~isempty(indj)K(2*i-1:2*i, 2*j-1:2*j) = K(2*i-1:2*i, 2*j-1:2*j) + eleK(2*indi-1:2*indi, 2*indj-1:2*indj, k);endendendendendFEMSolve.m文件：% 求解有限元方程function hFEM = FEMSolve(hFEM, vertdisind, R)K = hFEM.K;hFEM.vertdis = SolveFEM(vertdisind, R, K);for i = 1:hFEM.elecounthFEM.eletruct(i).vertdis(1:2:5) = hFEM.vertdis(2*hFEM.eleindex(i, :) - 1); % uhFEM.eletruct(i).vertdis(2:2:6) = hFEM.vertdis(2*hFEM.eleindex(i, :)); % vendendfunction vertdis = SolveFEM(vertdisind, R, K)vertdis = vertdisind;ind = find(vertdisind ~= 0);Krule = K(ind, ind); Rrule = R(ind);vertdis(ind) = Krule\Rrule; end。

MATLAB有限元分析与应用精选全文完整版

MATLAB有限元分析与应用精选全文完整版
function y = SpringElementForces(k,u)
%SpringElementForces This function returns the element nodal force
%
vector given the element stiffness matrix k
%
and the element nodal displacement vector u.
2019/11/28
§2-1 弹簧元
u1=U(1:2); f1=SpringElementForces(k1,u1);
f1 = -15.0000 15.0000
u2=U(2:3); f2=SpringElementForces(k2,u2);
f2 = -15.0000 15.0000
12
§3-1 弹簧元
%
modulus of elasticity E, cross-sectional
%
area A, and length L. The size of the
%
element stiffness matrix is 2 x 2.
y = [E*A/L -E*A/L ; -E*A/L E*A/L];
2019/11/28
3.1 单元刚度矩阵的形成
function y = SpringElementStiffness(k)
%SpringElementStiffness This function returns the element stiffness %matrix for a spring with stiffness k. %The size of the element stiffness matrix is 2 x 2.

在word中运行Matlab程序

在word中运行Matlab程序

在word中运行Matlab【1】在Matlab目录下，找到m-book.dot文件，双击运行【2】如果禁用宏，改成运行宏【3】这时打开的word，有一个notebook的菜单【4】选中要运行的程序【5】在notebook中，点击define input cell，这时程序变绿色了【6】再在notebook中，点击evaluate calc zone，就运行了。

///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// //////////////////////////使用matlab的notebook技术可以方便的实现这个功能，但是使用notebook的前提是matlab必须处在打开状态，因此主要用来制作科技文档。

装完Word和matlab后，在matlab主环境下运行“notebook –setup”看到“notebook setup is complete”就是安装成功了。

在Word的normal.dot模板相同目录下会有一个m-book.dot的文件，所谓在word中使用matlab其实就是加载这个模板，可以手工启动matlab，然后输入“notebook”，但这个方法比较麻烦。

可以在Word->工具->模板和加载项->模板中添加这个文件，这样word主菜单上就多了一个Notebook菜单，写一个命令后选中命令文字，然后Notebook->define input cell，指明这是一个输入，然后Notebook->Evaluate Cell就可以求解了。

与matlab中一样，命令后以分号结尾，不显示输出。

///////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////// ///////////////////////////////Matlab的Notebook软件工具设置及程序运行博战捷摘要介绍了在Matlab中如何设置Notebook软件工具，将Matlab程序嵌入中文Word。

(完整版)有限元大作业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程序

计算力学(有限元方法部分) 程序设计程序说明书程序1：平面问题的有限元分析文件名：h01.m算例文本：h01.txt输出文本：result1.txt使用前请先将h01.txt放入默认的文本读取路径（我的要求与m文件在同一文件夹内）！文本输入顺序：材料信息（编号、弹性模量、泊松比）注意：材料编号须按1、2、3、4……的顺序排列节点信息（编号、x坐标、y坐标）注意：节点编号须按1、2、3、4……的顺序排列约束信息（约束节点号、x方向有无约束、y方向有无约束、x方向位移、y 方向位移）有约束处填一正数，无约束处填0。

无约束处请勿输入位移。

单元信息（厚度、材料编号、节点编号，若为3节点单元，则第四个编号填0）注意：单元编号须按1、2、3、4……的顺序排列集中力（作用节点号、x方向力、y方向力）线荷载（作用边上的两个节点号、x方向分布力、y方向分布力）面荷载（作用单元号、x方向分布力、y方向分布力）程序可调部分：4-6行中可以指定输出哪些图像（按顺序依次为节点、单元图像，x、y方向位移图像，xx、yy、xy方向应力图像），第7行中可以指定输入的.txt文本名称。

文本输出内容：结点位移信息（节点号、x方向位移、y方向位移）单元形心处的应变信息（单元号、x方向正应变、y方向正应变、xy方向工程切应变）单元形心处的应力信息（单元号、x方向正应力、y方向正应力、xy方向切应力）本程序附有三角形单元自动加密前处理部分h01auto.m，其算例文本：h01coarse.txt，输出文本：h01new.txt。

它可以适用于本题的要求，在已给定本题所需全部信息的情况下将已有的单元加密为三角形单元。

其输出文本可直接作为上面程序的输入文本。

h01.m h01.txt h01auto.m h01coarse.txt欢迎交流与提问！附上邮箱：***************。

祝力学学习顺利！。

matlab有限元法

matlab有限元法

matlab有限元法
有限元法是一种常用的工程数值计算方法，广泛应用于结构力学、流体力学、热传导等领域。

它通过将复杂的连续体分割成有限个简单的单元，利用单元之间的相互关系来近似描述整个问题的解。

在工程实践中，有限元法已经成为一种不可或缺的分析工具。

有限元法的基本步骤包括建立数学模型、离散化、确定边界条件、求解方程、后处理等。

首先，需要将实际工程问题转化为数学模型，确定问题的几何形状、材料特性和载荷条件。

然后，将问题离散化，即将结构分割成有限个简单的单元，并确定单元之间的连接关系。

接下来，需要确定边界条件，即给定结构的边界约束和外部载荷。

然后，通过求解离散化后的方程组，得到问题的数值解。

最后，进行后处理，分析和展示结果。

有限元法的优点在于能够处理复杂的几何形状和边界条件，可以模拟各种不同的物理现象，并且具有较高的精度和可靠性。

它能够帮助工程师更好地理解和设计结构，提高工程的可靠性和安全性。

然而，有限元法也存在一些局限性。

首先，离散化过程会引入一定的误差，尤其是在模型中存在较大的变形或应力集中的情况下。

其次，求解大规模的方程组需要较高的计算资源和时间。

此外，有限元法对材料的本构关系和边界条件的设定要求较高，需要进行合理的模型假设和参数选择。

总的来说，有限元法是一种强大而灵活的工程分析方法，能够帮助工程师解决各种复杂的工程问题。

通过合理的模型建立和边界条件设定，以及精确的计算和后处理，可以得到准确可靠的结果，为工程设计和优化提供有力支持。

利用MATLAB生成Word文档

利用MATLAB生成Word文档

利用MATLAB生成Word文档fun ctio n ceshi_Word%利用MATLAB生成Word文档% ceshi_Word%% Copyright 2009 - 2010 xiezhh.% $Revisio n: 1.0.0.0 $ $Date: 2009/12/14 12:37:00 $%设定测试Word文件名和路径测试.doc'];%判断Word是否已经打开，若已打开，就在打开的Word中进行操作，否则就打开try%若Word服务器已经打开，返回其句柄WordWord = actxGetR unnin gServer('Word.Applicati on');catch%创建一个Microsoft Word服务器，返回句柄WordWord = actxserver('Word.Applicatio n');en d;%设置Word属性为可见Word.Visible = 1; % 或set(Word, 'Visible', 1);%若测试文件存在，打开该测试文件，否则，新建一个文件，并保存，文件名为测试if exist(filespec_user,'file');Docume nt = Word.Docume nts.Ope n(filespec_user);% Docume nt = in voke(Word.Docume nts,'Ope n',filespec_user);elseDocume nt = Word.Docume nts.Add;% Docume nt = in voke(Word.Docume nts, 'Add');Docume nt.SaveAs(filespec_user);endContent = Document.Content; % 返回Content 接口句柄Selection = Word.Selection; % 返回Selection 接口句柄Paragraphformat = Selection.ParagraphFormat; % 返回ParagraphFormat 接口句柄%页面设置Docume nt.PageSetup.TopMarg in = 60;Docume nt.PageSetup.BottomMarg in = 45;Docume nt.PageSetup.LeftMarg in = 45;Docume nt.PageSetup.RightMarg in = 45;%设定文档内容的起始位置和标题Con te nt.Start = 0; %设置文档内容的起始位置Word .doc%上边距60磅%下边距45磅%左边距45磅%右边距45磅title ='试卷分析'；Conten t.Text = title;%输入文字内容 Conten t.F on t.Size =16 ; %设置字号为16 Co nten t.Fo nt.Bold =4 ; %字体加粗Con te nt.P aragraphs.Alig nment = 'wdAlig nParagraphCe nter ：Select ion. TypeParagraph; xueqi = ' ( 2009 — 2010 学年第一学期)’；Selection.Text = xueqi; % 在当前位置输入文字内容Selectio n.Fo nt.Size = 12; % 设置字号为 12 Selectio n.Fo nt.Bold = 0; % 字体不加粗Tables = Docume nt.Tables.Add(Selectio n.Ra nge,12,9); %返回第1个表格的句柄DTI = Document.T ables.ltem(1); % 或 DTI = Tables; %设置表格边框DTI.Borders.OutsideLi neStyle = 'wdLi neStyleSi ngle'; DTI.Borders.OutsideLi neWidth = 'wdL in eWidth150pt';DTI.Borders.I nsideLi neStyle = 'wdLi neStyleSi ngle'; DTI.Borders.I nsideLi neWidth = 'wdL in eWidth150pt'; DTI.Rows.Alig nment = 'wdAlig nRowCe nter';DTI.Rows.ltem(8).Borders.ltem(1).Li neStyle = 'wdL in eStyleNo ne'; DTI.Rows.ltem(8).Borders.lte m( 3).Li neStyle = 'wdL in eStyleNo ne';DTI.Rows.Item(11).Borders.Item(1).Li neStyle = 'wdLi neStyleNo ne'; DTI.Rows.ltem(11).Borders.ltem(3).Li neStyle = 'wdLi neStyleNo ne'; %设置表格列宽和行高column_width = [53.7736,85.1434,53.7736,35.0094,...35.0094,76.6981,55.1887,52.9245,54.9057];% 定义列宽向量row_height = [28.5849,28.5849,28.5849,28.5849,25.4717,25.4717,...32.8302,312.1698,17.8302,49.2453,14.1509,18.6792]; % 定义行高向量%通过循环设置表格每列的列宽 for i = 1:9DTI.Colu mn s.Item(i).Width = colu mn _width(i);%居中对齐Selectio n.Start = Conten t.e nd; %设定下面内容的起始位置%回车，另起一段 Selectio n.MoveDow n;%光标下移（取消选中）paragraphformat.Alig nment = 'wdAlig nParagraphCe nter ：Select ion. TypeParagraph; %回车，另起一段Select ion. TypeParagraph;%回车，另起一段Selection.Font.Size = 10.5; % 设置字号为 10.5%居中对齐%在光标所在位置插入一个12行9列的表格end%通过循环设置表格每行的行高for i = 1:12DTI.Rows.ltem(i).Height = row_height(i);end%通过循环设置每个单元格的垂直对齐方式for i = 1:12for j = 1:9DTI.Cell(i,j).VerticalAlig nment = 'wdCellAlig nV erticalCe nter'; endend%合并单元格DTI.Cell(1,4).Merge(DTI.Cell(1,5));DTI.Cell(2, 4).Merge(DTI.Cell(2, 5));DTI.Cell(3, 4).Merge(DTI.Cell(3, 5));DTI.Cell(4, 4).Merge(DTI.Cell(4, 5));DTI.Cell(5, 2).Merge(DTI.Cell(5, 5));DTI.Cell(5, 3).Merge(DTI.Cell(5, 6));DTI.Cell(6, 2).Merge(DTI.Cell(6, 5));DTI.Cell(6, 3).Merge(DTI.Cell(6, 6));DTI.Cell(5, 1).Merge(DTI.Cell(6, 1));DTI.Cell(7, 1).Merge(DTI.Cell(7, 9));DTI.Cell(8, 1).Merge(DTI.Cell(8, 9));DTI.Cell(9, 1).Merge(DTI.Cell(9, 3));DTI.Cell(9, 2).Merge(DTI.Cell(9, 3));DTI.Cell(9, 3).Merge(DTI.Cell(9, 4));DTI.Cell(9, 4).Merge(DTI.Cell(9, 5));DTI.Cell(10, 1).Merge(DTI.Cell(10, 9));DTI.Cell(11,5).Merge(DTI.Cell(11,9));DTI.Cell(12, 5).Merge(DTI.Cell(12, 9));DTI.Cell(11, 1).Merge(DTI.Cell(12, 4));Selection.Start = Content.end; % 设置光标位置在文档内容的结尾Select ion. TypeParagraph; %回车，另起一段Selection.Text ='主管院长签字：年月日'；%输入文字内容Paragraphformat.Alig nment = 'wdAlig nParagraphRight'; % 右对齐Selection.MoveDown; % 光标下移%写入表格内容DTI.Cell(1,1).Range.Text ='课程名称'；DTI.Cell(1,3).Range.Text ='课程号'；DTI.Cell(1,5).Range.Text ='任课教师学院’；DTI.Cell(1,7).Range.Text ='任课教师DTI.Cell(2,1).Range.Text ='授课班级DTI.Cell(2,3).Range.Text ='考试日期DTI.Cell(2,5).Range.Text ='应考人数DTI.Cell(2,7).Range.T ext ='实考人数DTI.Cell(3,1).Range.T ext ='出卷方式DTI.Cell(3,3).Range.Text ='阅卷方式'; DTI.Cell(3,5).Range.T ext ='选用试卷A/B ：DTI.Cell(3,7).Range.Text ='考试时间'; DTI.Cell(4,1).Range.Text ='考试方式'; DTI.Cell(4,3).Range.T ext ='平均分'；DTI.Cell(4,5).Range.T ext ='不及格人数DTI.Cell(4,7).Range.Text ='及格率'；DTI.Cell(5,1).Range.T ext ='成绩分布'; DTI.Cell(5,2).Range.Text = '90 分以上DTI.Cell(5,3).Ra nge.Text = '80---89 分DTI.Cell(6,2).Ra nge.Text = '70--79 分DTI.Cell(6,3).Ra nge.Text = '60---69 分DTI.Cell(7,1).Ra nge.Text =['试卷分析(含是否符合教学大纲、难度、知识覆'盖面、班级分数分布分析、学生答题存在的共性问题与知识掌握情况、教学中 '存在的问题及改进措施等内容) '];DTI.Cell(7,1).Ra nge.ParagraphFormat.Alig nment = 'wdAlig nParagraphLeft ：DTI.Cell(9,2).Range.Text ='签字：：DTI.Cell(9,4).Range.Text ='年月日'； DTI.Cell(10,1).Range.Text ='教研室审阅意见：’；DTI.Cell(10,1).Ra nge.ParagraphFormat.Alig nment = 'wdAlig nParagraphLeft'; DTI.Cell(10,1).VerticalAlig nment = 'wdCellAlig nVerticalTop'; DTI.Cell(11,2).Range.Text ='教研室主任(签字)：年月日';DTI.Cell(11,2).Ra nge.ParagraphFormat.Alig nment = 'wdAlig nParagraphLeft'; DTI.Cell(8,1).Ra nge.P aragraphFormat.Alig nment = 'wdAlig nParagraphLeft'; DTI.Cell(8,1).VerticalAlig nment = 'wdCellAlig nV erticalTop'; DTI.Cell(9,2).Borders.ltem(2).L in eStyle = 'wdLi neStyleN on e'; DTI.Cell(9,2).Borders.ltem(4).L in eStyle = 'wdLi neStyleN on e'; DTI.Cell(9,3).Borders.Item(4).L in eStyle = 'wdLi neStyleN on e'; DTI.Cell(11,1).Borders.Item(4).Li neStyle = 'wdL in eStyleNo ne'; %如果当前工作文档中有图形存在，通过循环将图形全部删除Shape = Document.Shapes; % 返回Shapes 接口的句柄ShapeCount = Shape.Count; % 返回文档中Shape 对象的个数if ShapeCou nt ~= 0;for i = 1:ShapeCo unt;Shape.Item(1).Delete; % 删除第 1 个 Shape 对象 en d;en d;%产生标准正态分布随机数，画直方图，并设置图形属性人占%'; 人占 %'; 人占%'; 人占%';zft = figure(' un its',' no rmalized','positi on',...[0.280469 0.553385 0.428906 0.251302],'visible','off); % 新建图形窗口，设为不可见set(gca,'position',[0.1 0.2 0.85 0.75]); % 设置坐标系的位置和大小data = normrnd(0,1,1000,1); % 产生标准正态分布随机数hist(data); %绘制标准正态分布随机数的频数直方图grid on; % 添加参考网格xlabel('考试成绩'); %为X轴加标签ylabel('人数'); %为Y轴加标签%将图形复制到粘贴板hgexport(zft, '-clipboard');%将图形粘贴到当前文档里(表格的第8行第1列的单元格里)，并设置图形版式为浮于文字上方% Select ion .Ra nge.PasteSpecial;DTI.Cell(8,1).Ra nge. Paragraphs.ltem(1).Ra nge. PasteSpecial;Shape.ltem(1).WrapFormat.Type = 3;Shape.ltem(1).ZOrder('msoBringlnFrontOfText'); % 设置图片叠放次序为浮于文字上方delete(zft); %删除图形句柄Document.ActiveWindow.ActivePane.View.Type = 'wdPrintView'; % 设置视图方式为页面Document.Save; % 保存文档口Love is n ot a maybe thing. You know whe n you love some one.。

MATLAB语言基础 ch4在Word环境下使用MATLAB

MATLAB语言基础 ch4在Word环境下使用MATLAB

西安邮电学院计算机系
Matlab程序设计基础 Matlab程序设计基础
3．删去M-book文档中所有输出细胞 ．删去 文档中所有输出细胞 Notebook菜单项中的 菜单项中的Purge Output Cells 菜单项中的 命令可以删去M-book文档中所有输出细胞。 文档中所有输出细胞。 命令可以删去 文档中所有输出细胞 4．细胞转化为文本 ． 细胞转化为文本的方法是：选定细胞， 细胞转化为文本的方法是：选定细胞， 再选择Notebook菜单中的 菜单中的Undefine Cells命令。 命令。 再选择 菜单中的 命令 或将光标置于细胞之中，按组合键Alt＋ 。 或将光标置于细胞之中，按组合键 ＋U。
西安邮电学院计算机系
Matlab程序设计基础 Matlab程序设计基础
4.2 细胞的使用
4.2.1 输入输出细胞 1．输入细胞 ． 定义输入细胞的方法是：首先选中所需命令， 定义输入细胞的方法是：首先选中所需命令，然 后在Notebook菜单项中选择 菜单项中选择Define Input Cell命令， 命令， 后在 菜单项中选择 命令 于是被选中的MATLAB命令成为输入细胞。定义 于是被选中的 命令成为输入细胞。 命令成为输入细胞 输入细胞也可以在选中所需命令后， 输入细胞也可以在选中所需命令后，直接按组合 键Alt+D。 。 为了执行输入细胞，应选择 为了执行输入细胞，应选择Notebook菜单项中的 菜单项中的 Evaluate Cell命令或直接按组合键 命令或直接按组合键Ctrl+Enter。 命令或直接按程序设计基础
4.1.2 Notebook的启动 的启动 启动Notebook有两种方法：从Word中启动或 有两种方法： 启动 有两种方法 中启动或 命令窗口启动。 从MATLAB命令窗口启动。 命令窗口启动 1．从MATLAB中启动 中启动Notebook ． 中启动 2．从Word中启动 中启动Notebook ． 中启动

有限元 含时薛定谔方程 matlab

有限元 含时薛定谔方程 matlab

有限元含时薛定谔方程matlab摘要：1.有限元方法概述2.含时薛定谔方程简介3.MATLAB 在解决含时薛定谔方程中的应用4.结论与展望正文：1.有限元方法概述有限元方法是一种数值计算方法，广泛应用于固体力学、流体力学、热传导等领域。

其基本思想是将待求解的连续体划分为有限个小的、简单的子区域，即有限元，通过求解每个子区域的局部方程，最后得到整个系统的解。

这种方法有效地避免了求解复杂的偏微分方程，将问题简化为求解一组线性或非线性代数方程，从而降低了求解难度。

2.含时薛定谔方程简介含时薛定谔方程是量子力学中描述粒子在时间演化过程中的波函数演化的基本方程，其形式为i(Ψ/t) = HΨ，其中i 为虚数单位，为约化普朗克常数，Ψ为波函数，t 为时间，H 为哈密顿算子。

含时薛定谔方程是量子力学领域中一个重要的研究课题，对于解决粒子在特定势场中的时间演化问题具有重要意义。

3.MATLAB 在解决含时薛定谔方程中的应用MATLAB 是一种功能强大的数学软件，可以进行矩阵运算、数据分析、可视化等多种操作。

在解决含时薛定谔方程问题中，MATLAB 可以辅助完成以下任务：(1) 编写程序，实现有限元方法求解含时薛定谔方程。

通过将含时薛定谔方程离散化为有限元形式，可以利用MATLAB 进行高效的矩阵运算，从而得到波函数的时间演化。

(2) 绘制波函数随时间的变化图像。

MATLAB 提供了丰富的绘图功能，可以直观地展示波函数随时间的变化规律，便于分析和观察。

(3) 分析波函数的性质。

MATLAB 可以方便地对波函数进行数值分析，例如计算波函数的模长、相位等性质，为研究粒子的动态行为提供理论依据。

4.结论与展望有限元方法和MATLAB 在解决含时薛定谔方程问题中发挥了重要作用。

通过将复杂的量子力学问题转化为数值计算问题，可以降低问题的求解难度。

同时，MATLAB 提供了丰富的工具和功能，使得求解过程更加高效和直观。

有限元的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生成Word文档

MATLAB生成Word文档

MATLAB⽣成Word⽂档⼀、创建Microsoft Word服务器1. 创建Microsoft Word服务器try% 若Word服务器已经打开，返回其句柄WordWord = actxGetRunningServer('Word.Application');catch% 创建⼀个Microsoft Word服务器，返回句柄WordWord = actxserver('Word.Application');end2. 设置对象属性设置Word服务器为可见状态set(Word, 'Visible', 1); %或Word.Visible = 1;⼆、建⽴Word⽂本⽂档1. 新建空⽩⽂档%调⽤Add⽅法建⽴⼀个空⽩⽂档，并返回其句柄DocumentDocument = Word.Documents.Add;2. 页⾯设置%查看PageSetup接⼝的所有属性Document.PageSetup.get%页⾯设置Document.PageSetup.TopMargin =60;%上边距60磅Document.PageSetup.BottomMargin =45;%下边距45磅Document.PageSetup.LeftMargin =45;%左边距45磅Document.PageSetup.RightMargin =45;%右边距45磅%查看枚举类型属性VerticalAlignment的属性值Document.PageSetup.set('VerticalAlignment')3. 写⼊⽂字Content接⼝Content接⼝有很多属性和⽅法，可通过Content.get和Content.methodsview命令查看；Start属性⽤来获取或设定⽂字内容的起始位置，End属性⽤来获取或设定⽂字内容的终⽌位置，Text属性⽤来写⼊⽂字内容，Font属性⽤于字体设置，Paragraphs属性⽤于段落设置。

有限元的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有限元计算有限元计算是一种常用的数值计算方法，广泛应用于工程领域。

而Matlab作为一种强大的数值计算软件，提供了丰富的工具和函数，使得有限元计算变得更加简单和高效。

有限元计算是一种将连续问题离散化为有限个简单子问题的方法。

它将复杂的连续问题转化为离散的有限元网格，然后通过求解每个单元上的方程，最终得到整个问题的解。

有限元计算可以用于求解结构力学、流体力学、热传导等各种物理问题。

在Matlab中进行有限元计算，首先需要构建有限元模型。

有限元模型由节点和单元组成，节点是问题的离散点，单元是连接节点的基本单元。

在Matlab中，可以使用函数如meshgrid、linspace等来生成节点坐标，使用函数如delaunay、trimesh等来生成单元。

然后，需要定义问题的边界条件和加载条件。

边界条件是指在问题的边界上给定的约束条件，加载条件是指在问题中施加的外部力或位移。

在Matlab中，可以使用函数如boundary、findEdges等来定义边界条件，使用函数如force、displacement等来定义加载条件。

接下来，需要定义问题的材料性质和单元特性。

材料性质是指问题中所使用的材料的力学性质，单元特性是指单元的几何形状和材料性质。

在Matlab中，可以使用函数如materialProperties、elementProperties等来定义材料性质和单元特性。

然后，需要建立有限元方程。

有限元方程是通过对每个单元上的方程进行组装得到的整体方程。

在Matlab中，可以使用函数如stiffnessMatrix、loadVector等来建立有限元方程。

最后，需要求解有限元方程。

在Matlab中，可以使用函数如solve、eigs等来求解有限元方程。

求解得到的结果可以用于分析问题的应力、位移、变形等。

除了上述基本步骤，Matlab还提供了丰富的后处理工具和函数，用于可视化和分析有限元计算的结果。

matlab 有限元法

matlab 有限元法

matlab 有限元法
Matlab中的有限元法（Finite Element Method，FEM）是一种常用的数值分析方法，用于模拟和解决包括结构力学、热传导、流体力学等问题。

它将连续介质划分为离散的有限单元，通过建立数学模型和使用近似解法来求解。

下面是一般步骤来使用Matlab进行有限元分析：
1. 剖分网格：将要模拟的连续介质划分为离散的有限单元（如三角形或四边形元素）。

2. 建立数学模型：根据具体问题的物理方程或导引方程，建立线性或非线性的方程模型。

3. 施加边界条件：确定并施加边界条件，如位移、载荷或约束等。

4. 组装刚度矩阵和载荷向量（Assembly）：通过元素刚度矩阵的组装，得到总系统的刚度矩阵和载荷向量。

5. 求解方程：通过求解总系统的线性方程组，得到未知位移或其他需要的结果。

6. 后处理结果：对求解结果进行可视化或分析，如绘制应力分布、位移云图、应变曲线等。

Matlab提供了丰富的工具箱和函数，用于各种结构和物理问题的有限元分析，例如Partial Differential Equation Toolbox（部分微分方程工具箱）和Structural Analysis T oolbox（结构分析工具箱），其中包含了常用的有限元分析函数和设置界面。

另外，Matlab还支持用户自定义编程，允许使用脚本或函
数来实现特定的有限元算法。

总之，通过Matlab的有限元分析工具和编程能力，可以方便地进行各种结构和物理问题的数值分析和模拟。

matlab有限元编程

matlab有限元编程

是的，MATLAB可以用于有限元编程，用于解决各种结构、固体力学、热传导、电磁场等物理问题的数值模拟。

有限元法是一种数值方法，用于求解偏微分方程，特别适用于复杂的几何形状和边界条件的问题。

以下是一些使用MATLAB进行有限元编程的基本步骤：
1. 建立几何模型：首先，需要定义模型的几何形状和边界条件。

这包括确定节点和单元（元素）的位置。

2. 网格生成：将模型划分为离散的节点和单元。

这些节点和单元构成了有限元网格，可以使用MATLAB中的函数或者专用的网格生成工具来实现。

3. 定义材料属性和载荷：为每个单元分配适当的材料属性，如弹性模量、密度等。

同时，还要定义在结构上施加的力或边界条件。

4. 组装刚度矩阵和载荷向量：根据有限元法原理，将每个单元的局部刚度矩阵组装成全局刚度矩阵，同时将载荷向量组装成全局载荷向量。

5. 施加边界条件：根据边界条件，在全局刚度矩阵和载荷向量中施加约束条件。

6. 求解方程：通过求解线性方程组，得到节点的位移和其他所需的结果。

7. 后处理：根据求解结果，进行结果的可视化和分析。

可以绘制应力、应变分布图，计算位移、反应力等。

MATLAB提供了丰富的工具箱和函数来进行有限元编程，如Partial Differential Equation Toolbox（偏微分方程工具箱）、Finite Element Analysis Toolbox（有限元分析工具箱）等，可以大大简化有限元编程的过程。

在使用这些工具时，可以参考MATLAB的官方文档和例子，以及相关的有限元理论和方法。

matlab有限元分析

matlab有限元分析

matlab有限元分析MATLAB有限元分析，使用数值分析技术来解决工程及物理系统中的复杂问题，是一项重要的数学分析方法。

随着科学技术的发展，MATLAB有限元分析已经成为当今工程领域中不可或缺的研究手段之一。

1.简介MATLAB有限元分析是一种通过使用数学工具（如椭圆面元、半球元、三角元）和计算机技术来求解物理运动系统的中心技术。

它主要应用于解决结构力学、流体力学、热传导、电磁学和等离子体问题，以及其他复杂的工程学问题。

MATLAB有限元分析是一种快速应用,可以解决各种复杂的物理运动系统问题，特别是在计算机模拟问题上表现极为出色。

2.MATLAB有限元分析的应用MATLAB有限元分析可以用来解决许多复杂的物理系统问题，如：结构力学难题、流动或传热问题、电磁学问题和等离子体问题。

有限元分析可以用来解决任何物理系统的复杂性问题，以及复杂结构的动力学分析。

（1）结构力学问题MATLAB有限元分析可以被用于结构力学问题的求解，如航空航天和电力系统中的复杂结构分析，桥梁和石油等的结构力学分析，以及车辆的结构力学分析。

结构力学问题的解决，通过使用有限元分析，所获得的数值结果可以有效地利用计算机软件解决。

（2）流动或传热问题MATLAB有限元分析可以用于流体力学和热传导问题的研究，如医疗器械中的流动分析，原子力发电中的流动分析，航天器中的热传导分析，航空航天系统的热管传热分析等。

（3）电磁学问题MATLAB有限元分析也可以应用于电磁学问题，广泛应用于电磁场分析、电磁测量、电磁吸引及其他物理学问题。

这些应用以及在电子学中的各种应用，如电子元件的模拟分析、参数测量以及运动控制等，都得益于MATLAB有限元分析技术。

（4）等离子体问题MATLAB有限元分析技术也可以应用于等离子体系统的研究，广泛用于地球物理学、化学反应器、燃料电池、太空探索等复杂问题的研究。

3.MATLAB有限元分析的优势MATLAB有限元分析的优势是，可以解决复杂的物理系统问题，特别是可以有效地解决计算机模拟问题。

有限元matlab

有限元matlab
Richtmyer多步格式算出的结果并不理 想，不但左边有波动，而且光滑性也不好。 拉克斯-温德洛夫多步格式算出的结果比较不 错，虽然左边有点小波动，但是初始函数的 宽度和高度都保持的不错。 MacCormack多 步格式求得的结果和拉克斯-温德洛夫多步格 式算出的结果差不多。
38
双曲线偏微分方程——二维对流方程
A(i,i+1)=1;
%与边界无关的点离散
A(i,i+nx-2)=1;
A(i,i-nx+2)=1;
end
end
end
end
end
u=A\b;
%求线性方程组的解u
format short
11
五点差分格式算例求解
在命令窗口输入程序：u=peEllip5(51,0,2,51,0,2);
12
五点差分格式算例求解结果
如果网格更密的话，即采用更多的节点进行计算，会得到 更光滑的曲面。
13
工字型差分格式
14
工字型差分格式
注意：这里给出的边界仍是矩形边界；并且做网格剖分时要保证x 方向和y方向上的网格步长相等。
15
工字型差分格式
给出同样的算例求解以下拉普拉斯方程：
计算结果比较：用工字型差分格式计算得到的结果与五点差分格式得 到的结果差不多，但是前者角点上的算法不太好，这是格式自身的缺 陷。
for j=1:(n+M) u0(j)=IniU(minx+(j-M-1)*h); %向左延拓M个节点的函数值
end else
for j=1:(n+M) u0(j)=IniU(minx+(j-1)*h); %向左延拓M个节点的函数值
end end u1=u0; for k=1:M
  1. 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

MATLAB报告
Matlab程序求解简要过程如下：
（1）求取单元节点位移提取矩阵T
单元节点位移提取矩阵T本质上是置换矩阵群中的一个，结果可将任意杂乱的节点顺序置换成统一的顺序。

另一方面其作用是对单元刚度矩阵进行“升维操作”，将单元刚度矩阵统筹到整体刚度矩阵上来，便于对总体节点位移矩阵和支座反力进行求取。

本程序分析过程中对单元1的节点提取是按顺序编号1-2-3，对单元2的节点提取是按顺序编号2-3-4。

单元1的节点位移提取矩阵如下：
单元2的节点位移提取矩阵如下：
（2）求取单元几何矩阵B
单元1的节点按编号顺序1-2-3分别进行对几何函数矩阵或算子矩阵的bi逆时针操作，对ci顺时针操作；单元2的节点按编号顺序2-3-4分别进行对几何函数矩阵的bi顺时针操作，对ci逆时针操作.在MATLAB程序中通过mod()取模函数来达到对节点的顺时针或逆时针循环操作。

单元1的几何矩阵如下：
单元2的几何矩阵如下：
（3）求取应力矩阵S
单元应力矩阵满足S=D*B,其中D为弹性矩阵，B为单元几何矩阵
各单元的弹性矩阵如下：
单元1的应力矩阵如下：
单元2的应力矩阵如下：
（4）求取单元刚度矩阵K
单元刚度矩阵K满足公式K=B’*D*B*t*A,其中t为平面板的厚度，A为单元面积，且单元刚度矩阵为对称矩阵。

单元1的刚度矩阵如下：
单元2的刚度矩阵如下：
（5）求取总体刚度矩阵sumKK
由上述步骤求得的单元刚度矩阵K利用单元虚功原理和刚度方程可导出K’*δ=f，其中δ为单元节点位移列阵，f为单元等效节点载荷列阵，为了能将各个单元刚度方程统一到一个整体，便需要步骤（1）的单元节点提取矩阵对单元刚度方程进行变换，将两个变换结果联立便得到总体刚度方程，其中也可得到总体刚度矩阵sumKK，且总体刚度矩阵可由sumKK=ΣT’*K*T求得。

总体刚度矩阵如下：
（6）求取总体节点位移矩阵和支座反力
利用上述步骤提到的总体刚度方程sumKK*delta=F，其中delta为总体节点位移矩阵，F为总体等效节点载荷列阵。

利用总体边界条件可将总体刚度方程分成两部分，一部分先求取总体节点位移矩阵delta，最后代入总体刚度方程求取总体等效节点载荷列阵F，由此便得到支座反力。

总体节点位移矩阵如下：
支座反力如：
（7）求取单元应力
各单元应力满足θ=S*δ，其中S为单元应力矩阵，δ为单元节点位移列阵。

单元1的应力：
单元2的应力：
下图为应变前后的各单元对比图：
（注：可编辑下载，若有不当之处，请指正，谢谢!）。

相关文档
最新文档