(完整版)有限元大作业matlab---课程设计例子
有限元法MATLAB作业
例题：如图所示的受力均匀分布载荷作用的薄平板结构，将平板离散成两个线性三角元，如图所示，假定2210,0.3,0.025,3000/m E GPa v t m w KN ====。
求 （1）该结构的整体刚度矩阵。
（2）节点2和节点3的水平位移和垂直位移。
（3）节点1和节点4的支反力。
（4）每个单元的应力。
（5）每个单元的主应力和主应力方向角。
图2. 用双线性三角元离散化的薄木板 解：（1）离散化域我们将平板分为两个单元，4个节点，如图2所示，分布载荷的总作用力平均分给节点2和节点3，由于结构是薄平板，所以假定其属于平面应力情况。
MATLAB 中采用的计算单位是KN 和m 。
表1给出了该题的单元连通性。
（2）单元刚度矩阵通过调用MATLAB 的LinearTriangleElementStiffness 函数，得到两个单元矩阵k 1和k 2，每个矩阵都是66⨯的。
k1=gangdujuzhen(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)的源程序：function k1= gangdujuzhen( f,E,NU,t,xi,yi,xj,yj,xm,ym,p ) %UNTITLED4 此处显示有关此函数的摘要% 此处显示详细说明k1=gangdujuzhen(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)图1. 薄平板结构k1= f(E,NU,t,xi,yi,xj,yj,xm,ym,p);>>k1=gangdujuzhen(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0.25, 0, 0.25, 1)k1 =1.0e+06 *2.0192 0 0 -1.0096 -2.0192 1.00960 5.7692 -0.8654 0 0.8654 -5.76920 -0.8654 1.4423 0 -1.4423 0.8654-1.0096 0 0 0.5048 1.0096 -0.5048-2.0192 0.8654 -1.4423 1.0096 3.4615 -1.87501.0096 -5.7692 0.8654 -0.5048 -1.8750 6.2740k2=gangdujuzhen2(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)的源程序：function k2= gangdujuzhen2( f,E,NU,t,xi,yi,xj,yj,xm,ym,p )%UNTITLED3 此处显示有关此函数的摘要% 此处显示详细说明k2=gangdujuzhen2(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)k2= f(E,NU,t,xi,yi,xj,yj,xm,ym,p);>>k2=gangdujuzhen2(@LinearTriangleElementStiffness, 210000000, 0.3, 0.025, 0, 0, 0.5, 0, 0.5, 0.25, 1)k2 =1.0e+06 *1.4423 0 -1.4423 0.8654 0 -0.86540 0.5048 1.0096 -0.5048 -1.0096 0-1.4423 1.0096 3.4615 -1.8750 -2.0192 0.86540.8654 -0.5048 -1.8750 6.2740 1.0096 -5.76920 -1.0096 -2.0192 1.0096 2.0192 0-0.8654 0 0.8654 -5.7692 0 5.7692（3）集成整体刚度矩阵⨯的，因此为了得到整体刚度矩阵K，我们由于结构有4个节点，所以刚度矩阵是88⨯的零矩阵。
有限元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有限元分析与应用精选全文完整版
%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.
Matlab 有限元法计算分析程序编写
6) M函数文件 与命令文件不同，函数文件从外界只能看到传给它的输入 量和送出来的计算结果，而内部运作是看不见的。它的特 点是 (1)从形式上看，与命令文件不同，函数文件的第一行总是以 “function”引导的“函数申明行”。 (2)从运行上看，与命令文件运行不同，每当函数文件运行， MATLAB就会专门为它开辟临时工作空间，所有中间变量 都存放在函数工作空间中，当执行完文件最后一条指令和 遇到return时，就结束该函数文件的执行，同时该临时函数 工作空间及其所有的中间变量立即被清除。 (3)对于函数文件中的变量，如果不作特别说明，默认为临时 局部变量，这些临时变量就存放在函数的临时工作空间中， 当函数结束时他们被立即清除。与之相对应的是全局变量， 他们是通过global指令进行特别申明，这些全局变量可被几 个不同的函数共享。 • 函数文件的编辑也可用MATLAB editor/debugger。
有限元法计算分析程序编写
结构参数输入，包括
1）节点坐标值 2）单元类型以及连接信息 3）各单元的弹性模量、截面积（厚度）等 4）荷载形式以及作用位置、作用方向、荷载值 5）约束条件 6）输出信息
m j
对节点和单元分别编号 每个节点的自由度根据 节点号计算得到
i
y
o
x
计算结构的刚度矩阵
对各单元作如下的计算 a)计算单元刚度矩阵 b)计算坐标转换矩阵（如果需要） c)作坐标转换计算（如果需要） d)按自由度顺序叠加到总刚度矩阵中
MATLAB的使用方法
1) 最简单的计算器使用法 求[12＋2×(7－4)]÷32的算术运算结果 （1）用键盘在MATLAB指令窗中输入一下内容 (12+2*(7-4))/3^2 （2）在上述表达式输入完成后，按【Enter】键，该指令被执行 （3）在指令执行后，MATLAB指令窗中将显示一下内容 ans = 2 [说明] 加 + 减 乘 * 除 / 或 \ (这两个符号对于数组有不同的含义) 幂 ^ “ans”是answer的缩写，其含义是运算答案，它是MATLAB的一个默 认变量
《有限元基础教程》_【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有限元分析操作基础
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编程
end
精品..
正放四角锥网架定义
if e==1
定义网架上层节点
hu=input('输入网架上层节点行数'); %定义网架上层节点的行数 lu=input('输入网架上层节点列数'); %定义网架上层节点的列数
dis_xu=input('输入网架上层节点列间距'); %定义网架上层的行间距 dis_yu=input('输入网架上层节点行间距'); %定义网架上层的列间距
Element=zeros(21,2); for i=1:2:7
Element(5/2*i-3/2,:)=[i,i+1]; Element(5/2*i-1/2,:)=[i,i+2]; Element(5/2*i+1/2,:)=[i,i+3]; end for i=2:2:8 Element(5*i/2-1,:)=[i,i+1]; Element(5*i/2,:)=[i,i+2]; end Element(21,:)=[9,10];
end
荷载及边界条件
P=input(‘定义节点荷载，按[node1 P1;node2 P2;...]输入’); %网架荷载输入
BC=input(‘定义边界约束，按[node1 Conx Cony Conz;node2 Conx Cony Conz);...]输 入,Con代表x、y、z方向约束，取0为约束，取1无约束’); %网架边界条件
精品..
单元属性相同
if Cont1==0 AE1=input('请输入统一的截面面积与弹性模量，按[A E]输入');
AE=zeros(Msum,3);
matlab有限元分析实例
1.物理现象：这个对工程师来说是直观的物理现象和物理量，温度多少度，载荷是多大等等。
通常来说，用户界面中呈现的、用户对工程问题进行设置时输入的都是此类信息。
2.数学方程：将物理现象翻译成相应的数学方程，例如流体对应的是NS方程，传热对应的是传热方程等等；大部分描述这些现象的方程在空间上都是偏微分方程，偶尔也有ODE（如粒子轨迹、化学反应等）。
在这个层面，软件把物理现象“翻译”为以解析式表示的数学模型。
3.数值模型：在定义了数学模型，并执行了网格剖分后，商业软件会将数学模型离散化，利用有限元方法、边界元法、有限差分法、不连续伽辽金法等方法生成数值模型。
软件会组装并计算方程组雅可比矩阵，并利用求解器求解方程组。
这个层面的计算通常是隐藏在后台的，用户只能通过一些求解器的参数来干预求解。
有限元是一种数值求解偏微分方程的方法。
基本过程大致是设置形函数，离散，形成求解矩阵，数值解矩阵，后处理之类的。
MATLAB要把这些过程均自己实现，不过在数值求解矩阵时可以调用已有函数。
可以理解为MATLAB是一个通用的计算器，当然它的功能远不止如此。
而ANSYS之类的叫做通用有限元软件，针对不同行业已经将上述过程封装，前后处理也比较漂亮，甚至不太了解有限元理论的人也能算些简单的东西，当然结果可靠性又另说了。
比较两者，ANSYS之类的用起来容易得多，但灵活性不如MATLAB。
MATLAB用起来很困难，也有人做了一些模块，但大多数只能解决一些相对简单的问题。
对于大多数工程问题，以及某些领域的物理问题，一般都用通用有限元软件，这些软件还能添加一些函数块，用以解决一些需要额外设置的东西。
但是对于非常特殊的问题，以及一般性方程的有限元解，那只能用MATLAB或C，Fortran之类的了。
matlab课程设计大作业
matlab课程设计大作业一、教学目标本课程的教学目标是使学生掌握MATLAB基本语法、编程技巧以及MATLAB 在工程计算和数据分析中的应用。
通过本课程的学习，学生将能够熟练使用MATLAB进行简单数学计算、线性方程组求解、函数图像绘制等。
1.掌握MATLAB基本语法和编程结构。
2.了解MATLAB在工程计算和数据分析中的应用。
3.熟悉MATLAB的函数库和工具箱。
4.能够使用MATLAB进行简单数学计算。
5.能够使用MATLAB求解线性方程组。
6.能够使用MATLAB绘制函数图像。
7.能够利用MATLAB进行数据分析和处理。
情感态度价值观目标：1.培养学生对计算机辅助设计的兴趣和认识。
2.培养学生团队合作和自主学习的能力。
二、教学内容本课程的教学内容主要包括MATLAB基本语法、编程技巧以及MATLAB在工程计算和数据分析中的应用。
1.MATLAB基本语法：介绍MATLAB的工作环境、基本数据类型、运算符、编程结构等。
2.MATLAB编程技巧：讲解MATLAB的函数调用、脚本编写、函数文件编写等编程技巧。
3.MATLAB在工程计算中的应用：介绍MATLAB在数值计算、线性方程组求解、图像处理等方面的应用。
4.MATLAB在数据分析中的应用：讲解MATLAB在数据采集、数据分析、数据可视化等方面的应用。
三、教学方法本课程采用讲授法、案例分析法、实验法等多种教学方法相结合的方式进行教学。
1.讲授法：通过讲解MATLAB的基本语法、编程技巧以及应用案例，使学生掌握MATLAB的基本知识和技能。
2.案例分析法：通过分析实际工程案例，使学生了解MATLAB在工程计算和数据分析中的应用。
3.实验法：安排上机实验，使学生在实际操作中巩固所学知识，提高实际编程能力。
四、教学资源本课程的教学资源包括教材、实验设备、多媒体资料等。
1.教材：选用《MATLAB教程》作为主要教材，辅助以相关参考书籍。
2.实验设备：为学生提供计算机实验室，配备有MATLAB软件的计算机。
有限元课程设计matlab
有限元课程设计matlab一、课程目标知识目标：1. 学生能理解有限元分析的基本原理，掌握运用MATLAB进行有限元建模和求解的基本步骤。
2. 学生能够运用MATLAB软件进行简单物理场的有限元模拟，并解释模拟结果。
3. 学生掌握如何将实际问题抽象为有限元模型，并能够运用MATLAB进行模型参数的设定和调整。
技能目标：1. 学生能够独立操作MATLAB软件，进行有限元模型的构建和求解。
2. 学生能够通过MATLAB编程实现有限元模型的自动化处理，包括前处理、求解和后处理。
3. 学生通过解决实际问题，提高数值分析能力和计算机应用能力。
情感态度价值观目标：1. 学生培养对科学研究的兴趣，特别是在工程计算和仿真领域。
2. 学生通过解决实际问题，体会数学和工程结合的美，增强对工程问题的探究欲望。
3. 学生通过团队合作解决问题，培养协作精神和解决问题的能力。
本课程针对高年级本科生或研究生，他们具备一定的数学基础和编程能力。
课程性质偏重实践，旨在通过MATLAB这一工具将有限元理论应用于具体问题的求解。
课程目标旨在使学生不仅掌握理论知识，而且能够实际操作，将理论知识转化为解决实际问题的技能。
通过课程学习，学生应能够将所学知识应用于未来的学术研究或工程实践中。
二、教学内容1. 有限元方法基本原理回顾：包括有限元离散化、单元划分、形函数、刚度矩阵和载荷向量等概念。
- 教材章节：第二章 有限元方法基础2. MATLAB编程基础：介绍MATLAB的基本操作、数据结构、流程控制、函数编写等。
- 教材章节：第三章 MATLAB编程基础3. MATLAB中的有限元工具箱使用：学习如何使用MATLAB内置的有限元工具箱进行建模和求解。
- 教材章节：第四章 MATLAB有限元工具箱介绍4. 有限元模型构建与求解：结合实际问题，学习如何构建有限元模型，并进行求解。
- 教材章节：第五章 有限元模型构建与求解5. 实例分析与上机操作：通过案例分析，让学生实际操作MATLAB软件，解决具体的有限元问题。
计算力学(有限元)matlab编程大作业(空间网架)
逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：1 逐一输入第 i 个节点和其他节点的连接情况,若连接,输入 1,否则输入 0,令和自身连接情况为 0：0 输入结构约束总数：11 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：1 输入约束节点对应的位移号（输入范围：1~3） ） ：3 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：3 输入约束节点对应的位移号（输入范围：1~3） ） ：3 按照编号逐一输入结构约束编号（先输入节点号） ） ：4 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：4 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：1 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：2 按照编号逐一输入结构约束编号（先输入节点号） ） ：5 输入约束节点对应的位移号（输入范围：1~3） ） ：3 输入施加在结构上的总的力的个数：3 按照编号逐一输入力作用节点编号（先输入节点号） ） ：2 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：-1 输入力的大小：1000 按照编号逐一输入力作用节点编号（先输入节点号） ） ：2 逐一输入力对应的向量中的坐标值：0 逐一输入力对应的向量中的坐标值：1
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μμμμμμμ 题目：如图所示是一平面梁。
有限元方法与MATLAB程序设计 有限元分析可视化
2
A1 MATLAB图形输出函数
（2）空间图形
line(X,Y,Z) 画3维折线函数。X、Y和Z是3个一维数组，分别表示折线各结点坐标 plot3(X,Y,Z) 连接(X,Y,Z)坐标点的空间折线 patch(X,Y.Z,C) (X,Y,Z)为顶点的空间多边形，C为颜色 mesh(X,Y,Zxis equal;
% 设置轴比例相等
A2 桁架和刚架结构变形
s = 2;
if nargin>4, s = 3; end % 桁架为2；刚架为3
4
8
gm = max(abs(max(gxy)-min(gxy))); U = 2e-2*U*gm/max(abs(U)); G = gxy+[U(1:s:end),U(2:s:end)];
mesh(Z) 以Z矩阵列行下标为x,y轴自变量，用网线表示的曲面 surf(X,Y,Z) (X,Y,Z)坐标点张成的曲面
fsurf 以函数f(x,y)或x = funx(u,v), y = funy(u,v), z = funz(u,v)为参数的空间曲面 surf(Z) 以Z矩阵列行下标为x,y轴自变量画曲面 colorbar 图形中颜色对应的值
[1 1 1] w
白
[0 0 0] k
黑
表A2 线型设定符
propertyvalue '–' '--' ':' '–.'
'none'
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：MATLAB是美国MathWorks公司出品的商业数学软件，用于数据分析、无线通信、深度学习、图像处理与计算机视觉、信号处理、量化金融与风险管理、机器人，控制系统等领域。
MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室），软件主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。
它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式。
MATLAB和Mathematica、Maple并称为三大数学软件。
它在数学类科技应用软件中在数值计算方面首屈一指。
MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等。
MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。
在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。
重要功能·MATLAB®: MATLAB 语言的单元测试框架·Trading Toolbox™: 一款用于访问价格并将订单发送到交易系统的新产品·Financial Instruments Toolbox™: 赫尔-怀特、线性高斯和LIBOR 市场模型的校准和Monte Carlo 仿真·Image Processing Toolbox™: 使用有效轮廓进行图像分割、对10 个函数实现C 代码生成，对11 个函数使用GPU 加速·Image Acquisition Toolbox™: 提供了用于采集图像、深度图和框架数据的Ki nect® for Windows®传感器支持·Statistics Toolbox™: 用于二进制分类的支持向量机(SVM)、用于缺失数据的PCA 算法和Anderson-Darling 拟合优度检验·Data Acquisition Toolbox™: 为Digilent Analog Discovery Design Kit 提供了支持包·Vehicle Network Toolbox™: 为访问CAN 总线上的ECU 提供XCP。
有限元平面矩形单元MATLAB程序设计【范本模板】
有限元平面矩形单元MATLAB程序设计摘要本论文主要研究内容是有限元平面矩形单元的基本原理和MATLAB软件的图形用户界面及函数编程的基本知识，并根据有限元平面矩形单元的计算单元刚度矩阵,结点位移等原理和方法，运用MATLAB 语言编写弹性力学平面矩形单元的计算程序，使得程序能够完成对不同荷载类型，不同结构单元的计算，同时设计了图形用户界面(GUI),使得程序能够在图形用户界面（GUI)上对应的输入框内输入相应数据,并能够计算出单元刚度矩阵，等效结点荷载，有限元上结点位移和单元应力;最后将集中荷载和均布荷载情况下的实例进行对MATLAB计算得出的位移和单元结点应力与ANSYS计算值比较分析，验证其正确性和通用性,并总结了MATLAB软件在有限元的运用.关键词MATLAB程序设计;有限元；矩形单元；刚度矩阵;单元应力The Rectangle Unit Plane of Finite Element of theMATLAB ProgrammingAbstractThis dissertation mainly studies the radical principle of the rectangle unit plane of the finite element and the basic knowledge of the Graphical User Interface and the function programming of the MATLAB software 。
And make use of MATLAB language to make up the GUI and the computation program of MATLAB according to the theory and method of the rectangle unit plane of the finite element to compute the element stiffness matrix , nodes’displacement and so on about different types of loads and different elements 。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
有限元大作业程序设计学校：天津大学院系：建筑工程与力学学院专业：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[结论与体会]:通过本次的课程设计，我对有限元的概念有了更加深刻的理解，同时也弥补了平时学习是疏忽的地方，充实了有限元知识。