三角形常应变单元matlab程序的编制与使用

三角形单元刚度矩阵是指在有限元分析中用于描述三角形单元的刚度矩阵。

在有限元分析中，三角形单元是一种常用的离散单元，用于对实际结构进行离散化。

利用三角形单元刚度矩阵，可以对结构进行更精确的力学分析和计算。

1. 三角形单元简介三角形单元是有限元分析中最简单的单元之一。

它由三个节点组成，并且三角形内部的任意一点都可以由这三个节点线性插值得到。

三角形单元在有限元分析中应用广泛，尤其对于不规则结构的离散化具有很好的适用性。

2. 刚度矩阵概述在有限元分析中，结构的刚度矩阵描述了结构对外力的响应。

对于三角形单元而言，其刚度矩阵描述了该单元内部的节点之间的相互作用关系。

刚度矩阵的计算是有限元分析中的关键步骤，它直接影响了分析结果的精度和准确性。

3. 三角形单元刚度矩阵的计算三角形单元刚度矩阵的计算涉及到对单元内部的材料性质、几何形状和边界条件的考虑。

在MATLAB中，可以利用数值计算的工具和函数来进行三角形单元刚度矩阵的计算。

在进行计算时，需要考虑单元的形状函数、单元的形状函数导数和积分点的选择等因素，以确保计算结果的准确性。

4. MATLAB中的三角形单元刚度矩阵计算方法在MATLAB中，可以利用函数库或自定义函数来实现三角形单元刚度矩阵的计算。

通过编写相应的程序，可以根据单元的几何形状和材料性质计算出相应的刚度矩阵。

在计算过程中，需要考虑到MATLAB中的矩阵运算、数值积分和误差控制等技术，以确保计算结果的准确性和稳定性。

5. 刚度矩阵的应用三角形单元刚度矩阵的计算结果可以用于结构的稳定性分析、应力分析、变形分析等方面。

将刚度矩阵与外部载荷相乘，即可求解结构的位移响应。

通过对刚度矩阵的正交分解，还可以得到结构的固有频率和模态形态。

刚度矩阵在工程实践中具有重要的意义。

6. 结语三角形单元刚度矩阵是有限元分析中的重要概念，对于工程结构的分析和计算具有重要的意义。

在MATLAB中，通过编写相应的程序，可以实现对三角形单元刚度矩阵的准确计算。

matlab构造差分三角形（原创版）目录1.差分三角形的定义与作用2.MATLAB 中构造差分三角形的方法3.差分三角形的性质与应用正文1.差分三角形的定义与作用差分三角形（Difference Triangle）是一种用于数值计算的结构，主要用于离散差分和求和运算。

在实际应用中，差分三角形可用于数值积分、数值微分、线性方程组求解等。

差分三角形的核心思想是将一个函数在某一区间上的值用一组差分来表示，从而简化问题并降低计算复杂度。

2.MATLAB 中构造差分三角形的方法在 MATLAB 中，我们可以通过以下步骤构造差分三角形：（1）确定计算区间：首先，我们需要确定差分三角形的计算区间，通常用两个端点表示，例如：a 和 b。

（2）选择差分阶数：差分阶数决定了差分三角形的精度，通常选择为奇数。

常见的差分阶数有 1、3、5 等。

（3）编写差分矩阵：根据差分阶数，我们可以编写一个差分矩阵。

以奇数阶差分为例，假设差分阶数为 3，我们可以编写如下的差分矩阵：```matlab% 差分矩阵D = [1 -1 1 -1 1];```（4）计算差分三角形：根据给定的函数和差分矩阵，我们可以计算差分三角形。

以一个简单的函数为例：```matlab% 函数 f(x) = x^2f = @(x) x.^2;```我们可以通过以下代码计算差分三角形：```matlab% 计算差分三角形T = trapezoidal(f, a, b, D);```3.差分三角形的性质与应用差分三角形具有以下性质：（1）线性性质：对于任意两个函数 f(x) 和 g(x)，它们的差分三角形之和等于各自差分三角形之和，即 T(f+g) = T(f) + T(g)。

（2）常数性质：对于任意常数 c，它的差分三角形为 c*T(1)。

（3）积分性质：差分三角形可以用于数值积分。

例如，对于函数 f(x)，我们可以通过计算差分三角形的和来估计其定积分值。

（4）求解线性方程组：差分三角形还可以用于求解线性方程组。

实验三：三角形单元的形函数性质验证
一、 实验目的
1、加深对平面三角形单元有限元分析过程的理解；
2、掌握平面三角形单元形函数矩阵的求解过程和性质。

二、 实验要求
1、明确形函数矩阵的含义和求法；
2、根据题目要求，给出具体的计算过程；
3、编制相应的matlab 计算程序并调试运行。

三、 实验内容
用有限元法求图示平面三角形单元的形函数。

已知节点i,j,m 在xoy 平面中。

用所编程序验证形函数特性：
1、在单元任一点上，三个形函数之和等于1；
2、形函数N i 在i 点的函数值为1，在j 点及m 点的函数值为零；
3、三边上任一点的形函数与第三个顶点的坐标无关。

四、 实验提示
1、()() 21,y c x b a y x N i i i i ++∆=
，其中下标i , j , m 轮换; 2、()()m i j m i j i m m j j i y x y x y x y x y x y x ++++=∆2
1 -21; 3、j m i m m i j m m j i x x c y y b y x y x a -=-=-=;;，其中下标i , j , m 轮换。

matlab plot 填充的三角形我们需要明确绘制填充的三角形需要哪些元素：顶点坐标。

假设我们选择以下三个顶点作为示例：A(1,1)，B(3,1)，C(2,3)。

我们需要在MATLAB中定义这三个顶点的坐标。

可以使用以下代码来定义这三个点的坐标：A = [1,1];B = [3,1];C = [2,3];接下来，我们可以使用MATLAB的`fill`函数来绘制填充的三角形。

`fill`函数需要接受一个顶点坐标的矩阵作为输入，并使用这些坐标绘制一个填充的图形。

以下是使用`fill`函数绘制填充的三角形的代码：fill([A(1),B(1),C(1)],[A(2),B(2),C(2)],'b');在上述代码中，`[A(1),B(1),C(1)]`表示三个顶点的x坐标，`[A(2),B(2),C(2)]`表示三个顶点的y坐标，`'b'`表示填充颜色为蓝色。

我们可以将上述代码放在一个MATLAB脚本文件中，并运行该脚本文件来绘制填充的三角形。

绘制结果如下图所示：在绘制的图形中，我们可以看到一个填充的三角形，顶点分别为A(1,1)，B(3,1)，C(2,3)。

填充颜色为蓝色。

除了填充颜色，我们还可以使用其他颜色来填充三角形。

MATLAB提供了很多内置的颜色名称，例如红色（'r'），绿色（'g'），黄色（'y'），等等。

我们可以根据需要选择不同的颜色来填充三角形。

我们还可以使用MATLAB的`hold on`函数来保持当前图形，并在同一图形上绘制多个填充的三角形。

以下是一个示例代码，演示如何绘制多个填充的三角形：A1 = [1,1];B1 = [3,1];C1 = [2,3];A2 = [4,4];B2 = [6,4];C2 = [5,6];fill([A1(1),B1(1),C1(1)],[A1(2),B1(2),C1(2)],'r');hold on;fill([A2(1),B2(1),C2(1)],[A2(2),B2(2),C2(2)],'g');在上述代码中，我们首先定义了两个不同的三角形的顶点坐标（A1、B1、C1和A2、B2、C2），然后分别使用`fill`函数绘制这两个三角形，并分别使用红色和绿色来填充。

abaqus三角形单元创建
带有3个节点的三角形是最基础的二维有限元单元。

相比起其他单元，3节点三角形最大的特点是应变和应力在单元上的分布是相同的，所以也被称为Constant Strain Triangular (CST) 单元。

本文介绍如何通过用户自定义单元子程序UEL在ABAQUS环境中开发CST单元。

单元为平面应变单元，对标ABAQUS中的CPE3。

由于ABAQUS不支持显示用户自己开发的单元，所以需要通过UMAT辅助显示UEL计算得到的结果。

大致的流程为：
（1）在input中定义两种单元，一种是用户自己开发的“真”单元（UEL），另一种是“假”单元（abaqus内置的单元，但是材料用UMAT），用于后处理显示。

这两种单元享用相同的节点标号，但是有不同的单元标号；
（2）在UEL和UMAT中都定义COMMON block，比如: common/KUSER/ UVAR；
（3）UEL计算得到结果，将需要显示的结果赋值给SVARS （UEL中的场变量）；
（4）在UEL中把SVARS 的数值赋予UVAR；
（5）UMAT中的SDV接收UVAR，另外需要写一个弹性的本构，弹性模量 E 可以设置得很小.。

傻瓜运行程序，全自动过程首先建立一个数据输入的M文件：syms o;p=1;E=210e6;%输入弹性模量；NU=0.3;%输入泊松比；t=0.025;%输入材料厚度；x=[0,0.5;0.25,0.5;0.125,0.375;0,0.25;0.25,0.25;0.125,0.125;0,0;0.25,0;0.375,0.125;0.5,0.25;0.5,0]; %输入每个结点的坐标值；ff=[1,3,2;1,4,3;3,5,2;3,4,5;4,6,5;4,7,6;5,6,8;6,7,8;5,8,9;5,9,10;8,11,9;9,11,10];%输入每个单元的信息；F=[o;o;0;0;0;0;o;o;0;-12.5;0;0;o;o;0;0;0;0;0;0;0;0];%输入总体结点载荷矢量，其中未知量用符号变量'o'表示；U=[0;0;o;o;o;o;0;0;o;o;o;o;0;0;o;o;o;o;o;o;o;o];%输入总体结点位移矢量，其中未知量用符号变量'o'表示；然后建立一个运行程序的M文件：K=zeros(2*size(x),2*size(x));%生成初始整体刚度矩阵；sigma=zeros(size(ff),3);%生成初始单元应力矩阵；s=zeros(size(ff),3);%生成初始单元主应力应力及主应力方向角矩阵；for i=1:size(ff)%i表示第i个单元；kk=LinearTriangleElementStiffness(E,NU,t,x(ff(i,1),1),x(ff(i,1),2),x(ff(i,2),1),x(ff(i,2),2),x(ff(i,3), 1),x(ff(i,3),2),p);%计算每个单元的单元刚度矩阵；K=LinearTriangleAssemble(K,kk,ff(i,1),ff(i,2),ff(i,3));%集成整体刚度矩阵；end%根据每个单元的单元刚度矩阵，循环集成整体刚度矩阵；k=K;for i=1:size(U)if U(i)==0F(i)=0;k(i,:)=0;k(:,i)=0;k(i,i)=1;elseendend%将未知的结点载荷及其对应的总体刚度的行和列定义为0，并定义对应的i行i列为1；U=k\F;F=K*U;%计算总体结点载荷矢量for i=1:size(ff)l=ff(i,1);m=ff(i,2);n=ff(i,3);uu=[U(2*l-1);U(2*l);U(2*m-1);U(2*m);U(2*n-1);U(2*n)];%建立单元结点位移矢量；ss=LinearTriangleElementStresses(E,NU,x(ff(i,1),1),x(ff(i,1),2),x(ff(i,2),1),x(ff(i,2),2),x(ff(i,3),1), x(ff(i,3),2),p,uu);%计算单元应力；sss=LinearTriangleElementPStresses(ss);%计算单元主应力和主应力方向角；sigma(i,1)=double(ss(1,1));sigma(i,2)=double(ss(2,1));sigma(i,3)=double(ss(3,1));%将每个单元应力赋值到单元应力矩阵sigma中；s(i,1)=double(sss(1,1));s(i,2)=double(sss(2,1));s(i,3)=double(sss(3,1));%将每个单元主应力和主应力方向角赋值到单元主应力应力及主应力方向角矩阵s中；endU=double(U);F=double(F);uu=double(uu);K%以下程序输出结点信息到‘节点信息.txt’中；fid=fopen('C:\Documents and Settings\LXY\桌面\节点信息.txt','wt');fprintf(fid,'\n%s\n','-------------------------------------------------------节点信息----------------------------------------------------------');fprintf(fid,'\n%s\n',' 节点号x坐标y坐标x位移y位移水平支反力垂直支反力');for i=1:size(x)fprintf(fid,'%10d%18.8f%18.8f%18.8f%18.8f%18.8f%18.8f\n',i,x(i,1),x(i,2),U(2*(i-1)+1),U(2*(i-1)+2),F(2*(i-1)+1),F(2*(i-1)+2));endfclose(fid);%以下程序输出单元信息到‘单元信息.txt’中；fid=fopen('C:\Documents and Settings\LXY\桌面\单元信息.txt','wt');fprintf(fid,'\n%s\n','---------------------------------------------------------------单元信息----------------------------------------------------------------------');fprintf(fid,'\n%s\n',' 单元号结点1 结点2 结点3 x方向拉应力y方向拉应力剪应力主应力方向角'); for i=1:size(ff)fprintf(fid,'%10d%18.8f%18.8f%18.8f%18.8f%18.8f%18.8f%18.8f\n',i,ff(i,1),ff(i,2),ff(i,3),sigma( i,1),sigma(i,2),sigma(i,3),s(i,3));endfclose(fid);。

matlab函数文件的编写和调用
Matlab函数文件的编写和调用步骤如下：
1. 打开Matlab软件，进入Matlab命令窗口；
2. 输入命令“edit 函数名”，创建新的函数文件，或者在文件浏览器中手动创建一个.m格式的函数文件；
3. 在函数文件中编写函数代码，包括函数名、输入参数、输出参数和函数体等；
4. 使用保存按钮将函数文件保存在Matlab当前的工作目录中；
5. 在Matlab命令窗口中，使用函数名加上输入参数来调用该函数，函数将返回输出参数的值。

例如，编写一个计算两数之和的函数sum：
1. 打开Matlab软件，进入Matlab命令窗口；
2. 输入命令“edit sum”，弹出编辑窗口；
3. 在编辑窗口中输入以下代码：
function c = sum(a, b)
% 计算两个数的和
c = a + b;
end
4. 保存函数文件“sum.m”在Matlab当前的工作目录中；
5. 在Matlab命令窗口中，输入“sum(2, 3)”调用函数，函数将返回5。

需要注意的是，函数文件名必须与函数名相同，并且保存在Matlab当前的工作目录中或其他Matlab可以搜索到的路径中。

调用函数时，输入参数的个数和类型必须与函数定义中的输入参数一致，否则将报错或产生意外结果。

matlab的德劳内三角形-回复MATLAB的德劳内三角形德劳内三角形是计算几何中一个重要的概念，用于解决许多实际问题，比如计算点云的凸壳、三角网格生成等。

在MATLAB中，我们可以使用内置的函数来计算德劳内三角形。

本文将以德劳内三角形为主题，分步介绍如何使用MATLAB进行计算。

第一步：创建点集首先，我们需要创建一个点集。

在MATLAB中，可以使用数组来表示点，每个点可以由两个或三个坐标值表示。

假设我们要创建一个包含10个随机点的点集，可以使用如下代码：matlabnumPoints = 10; 点的数量points = rand(numPoints, 2); 随机生成2维点坐标这样就创建了一个10x2的数组points，其中每一行表示一个点的坐标。

第二步：计算德劳内三角形接下来，我们可以使用MATLAB的delaunay函数来计算德劳内三角形。

该函数将点集作为输入，并返回一个包含三角形的连接关系矩阵，其中每一行表示一个三角形的三个顶点索引。

代码如下：matlabtri = delaunay(points(:, 1), points(:, 2));tri是一个包含三个列的矩阵，每一行表示一个三角形的三个顶点索引。

第三步：绘制德劳内三角形我们可以使用MATLAB的triplot函数来绘制德劳内三角形。

该函数将连接关系矩阵和点集作为输入，并绘制出所有的三角形。

代码如下：matlabtriplot(tri, points(:, 1), points(:, 2));上述代码会绘制出所有的三角形，以及对应的点集。

第四步：可视化结果为了更好地显示德劳内三角形，我们可以使用MATLAB的scatter函数来绘制点集，并使用不同的颜色来表示不同的三角形。

代码如下：matlabscatter(points(:, 1), points(:, 2), 50, 'filled');hold on;triplot(tri, points(:, 1), points(:, 2));hold off;以上代码首先使用scatter函数绘制点集，其中第三个参数50表示点的大小，'filled'表示填充点的颜色。

第三章平面问题的有限元法本章通过三角形常应变单元，介绍有限元法应用于弹性体应力分析的基本原理和方法：包括弹性体的离散化，单元特性的分析，刚度矩阵的建立，等效节点力的计算，解答的收敛性以及实施步骤和注意事项，同时对形函数的性质也作简要的叙述。

第一节三角形常应变单元一、结构离散化用有限元法分析弹性力学平面问题，第一步就是把原来的连续的弹性体离散化。

(a) (b)将整个结构（平板）划分成有限个三角形。

这样的三角形称为单元（三角形单元）。

三角形单元的顶点取为节点。

3节点三角形单元用边界节点间的直线段来近似板的曲线边界。

这些三角形在其节点处相互连接，组成一个单元集合体，以代替原来的弹性体。

注：1. 全部节点和全部单元一般由1开始按自然顺序编号。

2. 节点编码：总码-----------用于整体分析，如1,2，…，按自然顺序编号局部码--------用于单元分析，i，j，m 要求按逆时针方向的顺序进行编码每个单元的节点局部码i，j，m和节点总码有一一对应的关系3. 单元间不能有重叠4. 一个单元的任一顶点不许为另一单元任一边的内点5. 所有作用在单元上的载荷，包括集中载荷、表面载荷和体积力，都按虚功等效的原则移置到节点上，成为等效节点载荷。

二、位移模式1. 单元节点位移列阵iu设单元e 的节点号码为i ，j ，m 。

由弹性力学平面问题可知，单元内任意一点有两个位移分量u ，v ，记为{}Tf uv ⎡⎤⎣⎦=故每个节点也有两个位移分量，因此称节点自由度为2。

节点i 位移分量记为{}Ti ii u v δ⎡⎤⎣⎦= （i ，j ，m 轮换）则3个节点的位移分量用列阵表示为{}i ei i e j j j m m m u v u v u v δδδδ⎧⎫⎪⎪⎪⎪⎧⎫⎪⎪⎪⎪⎪⎪⎪⎪⎨⎬⎨⎬⎪⎪⎪⎪⎪⎪⎪⎪⎩⎭⎪⎪⎪⎪⎩⎭==（3-1）称为单元节点位移列阵（向量）。

单元自由度是6。

2. 位移模式结构承受载荷以后发生位移，但位移分布事先并不知道。

matlab三角函数写法Matlab是一款广泛应用于数学和工程领域的数学软件，它提供了丰富的三角函数库，使得我们可以轻松地编写和求解三角函数问题。

本篇文章将详细介绍Matlab中常用的三角函数，以及如何正确使用它们进行绘图和计算。

一、正弦函数与余弦函数在Matlab中，正弦函数和余弦函数的符号分别为sin和cos。

它们的输入参数为角度（弧度制），输出为该角度的正弦值或余弦值。

例如，以下代码将绘制一个以0到2π为范围的圆，并标出每个角度对应的正弦值和余弦值：```matlabtheta = 0:pi/50:2*pi; % 生成角度向量sin_val = sin(theta); % 计算正弦值cos_val = cos(theta); % 计算余弦值plot(theta, sin_val); % 绘制正弦函数图像hold on; % 保持当前图像，以便绘制其他图形plot(theta, cos_val); % 绘制余弦函数图像legend('sin(x)', 'cos(x)'); % 添加图例说明```二、正切函数与余切函数Matlab中，正切函数和余切函数的符号分别为tan和cot。

它们的输入参数也为角度（弧度制），输出为该角度的正切值或余切值。

需要注意的是，在Matlab中，角度的范围通常默认为度数制，但为了与角度为弧度制的计算相匹配，我们在编写代码时需要注意单位转换。

三、其他三角函数除了正弦、余弦和正切函数外，Matlab还提供了其他一些常用的三角函数，如正弦的平方函数sin^2、反正切函数atan等。

这些函数的符号分别为sin^2、acos等。

在使用这些函数时，需要注意输入参数的单位和范围，以确保计算结果的准确性。

四、绘图与计算在Matlab中，我们可以使用plot函数绘制三角函数的图像，使用相关函数进行数值计算和统计分析。

例如，以下代码将绘制一个以0到π为范围的三角函数的图像，并使用相关函数进行数值计算：```matlabtheta = linspace(0, pi, 100); % 生成角度向量sin_val = sin(theta); % 计算正弦值cos_val = cos(theta); % 计算余弦值atan_val = atan(theta); % 计算反正切值plot(theta, sin_val, 'r'); % 绘制正弦函数图像，颜色为红色hold on; % 保持当前图像，以便绘制其他图形plot(theta, cos_val, 'b'); % 绘制余弦函数图像，颜色为蓝色plot(theta, atan_val, 'g'); % 绘制反正切函数图像，颜色为绿色legend('sin(x)', 'cos(x)', 'atan(x)'); % 添加图例说明```以上就是Matlab中三角函数的常用写法及绘图方法。

《有限元法》实验报告专业班级力学（实验）1601 姓名田诗豪学号 10提交日期实验一（30分）一、实验内容编写一个计算平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的MATLAB 函数文件[B3,S3,K3] = ele_mat_tri3(xy3,mat)，其中：输入变量xy3为结点坐标数组，mat 为材料参数矩阵；输出变量B3为应变矩阵，S3为应力矩阵，K3为单元刚度矩阵。

（要求给出3个不同算例进行验证，并绘制出单元形状和结点号）二、程序代码通用函数function [B3,S3,K3] = ele_mat_tri3(xy3,mat)%生成平面3结点三角形单元的应变矩阵、应力矩阵和单元刚度矩阵的功能函数%*********变量说明****************%xy3------------------结点坐标数组%mat------------------材料参数矩阵（弹性模量，泊松比，壁厚）%B3-------------------应变矩阵%S3-------------------应力矩阵%K3-------------------单元刚度矩阵%*********************************xyh=[1,xy3(1,1),xy3(1,2);1,xy3(2,1),xy3(2,2);1,xy3(3,1),xy3(3,2)];A=*det(xyh);A=abs(A);D=mat(1)/(1-mat(2)^2)*[1,mat(2),0;mat(2),1,0;0,0,(1-mat(2))/2];b=zeros(1,3);c=zeros(1,3);%*********************************for i=1:3if i==1j=2;m=3;elseif i==2j=3;m=1;elsej=1;m=2;endb(i)=xy3(j,2)-xy3(m,2);c(i)=xy3(m,1)-xy3(j,1);end%*********************************B31=1/(2*A)*[b(1),0;0,c(1);c(1),b(1)];B32=1/(2*A)*[b(2),0;0,c(2);c(2),b(2)];B33=1/(2*A)*[b(3),0;0,c(3);c(3),b(3)];B3=[B31,B32,B33];%*********************************S3=D*B3;%*********************************K3=A*mat(3)*B3'*D*B3;主程序clear;clc;%*********输入结点坐标数组********xy3=[0,0;5,1;1,4];mat=[3e6,,];%****输入材料参数矩阵(弹性模量，泊松比，壁厚)****[B3,S3,K3]=ele_mat_tri3(xy3,mat)三、算例分析算例1：如图1所示三角形单元，结点坐标为1(0,0)，2(5,2)，3(1,4)，弹性模量为200GPa，泊松比为、厚度为。

三角形单元有限元程序设计%*******************************************************************% NNODE 单元节点数NPION 总结点数% NELEM 单元数NVFIX 受约束边界点数% FIXED 约束信息数组NFORCE 节点力数% FORCE 节点力数组COORD 结构节点坐标数组% LNODS 单元定义数组% YOUNG 弹性模量POISS 泊松比% THICK 厚度 B 单元应变矩阵(3*6)% D 单元弹性矩阵(3*3) S 单元应力矩阵(3*6)% A 单元面积ESTIF 单元刚度矩阵% ASTIF 总体刚度矩阵ASLOD 总体荷载向量% ASDISP 节点位移向量ELEDISP 单元节点位移向量% STRESS 单元应力FP1 数据文件指针%*******************************************************************format short e %设定输出类型clear %清除内存变量FP1=fopen('1-1.txt','rt'); %打开输入数据文件存放初始数据%读入控制数据NELEM=fscanf(FP1,'%d',1), %单元个数（单元编码总数）NPION=fscanf(FP1,'%d',1), %结点个数（结点编码总数）NVFIX=fscanf(FP1,'%d',1) %受约束边界点数NFORCE=fscanf(FP1,'%d',1), %结点荷载个数YOUNG=fscanf(FP1,'%e',1), %弹性模量POISS=fscanf(FP1,'%f',1), %泊松比THICK=fscanf(FP1,'%f',1) %厚度LNODS=fscanf(FP1,'%d',[3,NELEM])' %单元定义数组（单元结点号）%相应为单元结点号（编码）、按逆时针顺序输入COORD=fscanf(FP1,'%f',[2,NPION])' %结点坐标数组%坐标：x，y坐标（共NPOIN组）FORCE=fscanf(FP1,'%f',[3,NFORCE])' %结点力数组%(n,3) n:受力结点数目,(n,1):作用点,(n,2):x方向,(n,3):y方向FIXED=fscanf(FP1,'%d',[3,NVFIX])' %约束信息数组% (n,1):约束点(n,2)与(n,3)分别为约束点x方向和y方向的约束情况,受约束为1否则为0%*******************************************************************%生成单元刚度矩阵并组成总体刚度矩阵ASTIF=zeros(2*NPION,2*NPION); %生成特定大小总体刚度矩阵并置0%*******************************************************************for i=1:NELEM%生成弹性矩阵DD= [1 POISS 0;POISS 1 0;0 0 (1-POISS)/2]*YOUNG/(1-POISS^2)%*******************************************************************%计算当前单元的面积A=det([1 COORD(LNODS(i,1),1) COORD(LNODS(i,1),2);1 COORD(LNODS(i,2),1) COORD(LNODS(i,2),2);1 COORD(LNODS(i,3),1) COORD(LNODS(i,3),2)])/2%*******************************************************************%生成应变矩阵Bfor j=0:2b(j+1)= COORD(LNODS(i,(rem((j+1),3))+1),2)-COORD(LNODS(i,(rem((j+2),3))+1),2);c(j+1)=-COORD(LNODS(i,(rem((j+1),3))+1),1)+COORD(LNODS(i,(rem((j+2),3))+1),1); endB=[b(1) 0 b(2) 0 b(3) 0;...0 c(1) 0 c(2) 0 c(3);...c(1) b(1) c(2) b(2) c(3) b(3)]/(2*A);B1( :,:,i)=B;%*******************************************************************%求应力矩阵S=D*BS=D*B;ESTIF=B'*S*THICK*A; %求解单元刚度矩阵a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号for j=1:3for k=1:3ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)...=ASTIF((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+ESTIF(j*2-1:j*2,k*2-1:k*2);%根据节点编号对应关系将单元刚度分块叠加到总刚矩阵中endendend%**************************************************************************%将约束信息加入总体刚度矩阵（对角元素改一法）for i=1:NVFIXif FIXED(i,2)==1ASTIF(:,(FIXED(i,1)*2-1))=0; %一列为零ASTIF((FIXED(i,1)*2-1),:)=0; %一行为零ASTIF((FIXED(i,1)*2-1),(FIXED(i,1)*2-1))=1;%对角元素为1endif FIXED(i,3)==1ASTIF( :,FIXED(i,1)*2)=0; %一列为零ASTIF(FIXED(i,1)*2,:)=0; %一行为零ASTIF(FIXED(i,1)*2 ,FIXED(i,1)*2)=1; %对角元素为1 endend%************************************************************************** %生成荷载向量ASLOD(1:2*NPION)=0; %总体荷载向量置零for i=1:NFORCEASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);end%************************************************************************** %求解内力ASDISP=ASTIF\ASLOD' %计算节点位移向量ELEDISP(1:6)=0; %当前单元节点位移向量for i=1:NELEMfor j=1:3ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);%取出当前单元的节点位移向量endiSTRESS=D*B1(:, :, i)*ELEDISP' %求内力endfclose(FP1); %关闭数据文件程序应用举例 134212yxFF F=100NF=100N如上图所示，矩形平板，一端固支，一端受集中力，划分为两个三角形单元，共4个结点。