平面三角形单元Matlab程序

合集下载

算例分析

一、平面3节点三角形单元分析的算例如图所示为一矩形薄平板，在右端部受集中力F =10000N 作用，材料常数为：弹性模量7110E Pa =⨯、泊松比13μ=，板的厚度为0.1t m =，试按平面应力问题计算各个节点位移及支座反力。

解：（1）结构的离散化与编号对该结构进行离散，单元编号及节点编号如图4-20(b)所示，即有二个3节点三角形单元。

载荷F 按静力等效原则向节点1、节点2移置等效。

节点位移列阵：[]11223344Tq u v u v u v u v =节点外载列阵：00000022TF F F ⎡⎤=--⎢⎥⎣⎦约束的支反力列阵：33440000Tx y x y R R R R R ⎡⎤=⎣⎦其中33(,)x y R R 和44(,)x y R R 分别为节点3和节点4的两个方向的支反力。

（2）各个单元的描述当两个单元取图示中的局部编码(i ,j ,m )时，其单元刚度矩阵完全相同，即(1),(2)ii ij im ji jj jm mi mjmm k k k k k k k k k k ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦=（3）建立整体刚度方程按单元的位移自由度所对应的位置进行组装可以得到整体刚度矩阵，该组装过程可以写成(1)(2)k k k =+具体写出单元刚度矩阵的各个子块在总刚度矩阵中的对应位置如下代入整体刚度方程Kq =P 中，有（4）边界条件的处理及刚度方程求解该问题的位移边界条件为33440,0,0,0u v u v ====，将其代入上式中，划去已知节点位移对应的第5行至第8行（列），有由上式可求出节点位移如下[]1122[] 1.888.99 1.508.42TT F u v u v Et=--- （5）支反力的计算将所求得的节点位移式代入总刚度方程中，可求得支反力如下3112924()23233x Et R u v v F =--+=- 31129214()0.0732333y Et R u v u F =--+=-42292()2323x Et R u v F =--=422921() 1.073233y Et R u v F =--=二、MATLAB —平面3节点三角形单元分析的算例(Triangle2D3Node)解：（1）结构的离散化与编号将结构离散为二个3节点三角形单元，单元编号及节点编号如图4-20(b)所示。

matlab旋转三角形的实现。

如何使用MATLAB旋转三角形在MATLAB中，我们可以使用矩阵变换的方法来实现三角形的旋转。

本文将介绍如何使用MATLAB进行三角形的旋转操作。

我们需要定义三角形的顶点坐标。

假设三角形的顶点坐标为A(x1, y1)，B(x2, y2)，C(x3, y3)。

可以使用MATLAB的矩阵表示来表示三角形的顶点坐标，例如：A = [x1, y1];B = [x2, y2];C = [x3, y3];顶点坐标可以根据实际情况进行定义。

接下来，我们需要定义旋转角度。

假设旋转角度为θ度。

可以使用MATLAB的角度转弧度函数来将角度转换为弧度，例如：theta = deg2rad(θ);其中，θ为旋转角度，theta为转换后的弧度。

然后，我们可以使用MATLAB的矩阵变换函数来进行三角形的旋转。

MATLAB中提供了一个名为makehgtform的函数，可以用于生成旋转矩阵。

旋转矩阵可以通过旋转角度和旋转轴来定义。

在本例中，我们使用makehgtform函数来生成绕Z轴旋转的变换矩阵，代码如下：rotation_matrix = makehgtform('zrotate', theta);其中，rotation_matrix为旋转矩阵。

接下来，我们需要将三角形的顶点坐标转换为齐次坐标。

可以通过在顶点坐标后添加一个额外的维度，并赋值为1来实现。

代码如下：A_homo = [A, 1];B_homo = [B, 1];C_homo = [C, 1];其中，A_homo，B_homo和C_homo分别为转换后的齐次坐标。

然后，我们可以使用MATLAB的矩阵乘法来进行三角形的旋转变换。

代码如下：A_rotated = A_homo * rotation_matrix;B_rotated = B_homo * rotation_matrix;C_rotated = C_homo * rotation_matrix;其中，A_rotated，B_rotated和C_rotated为旋转后的顶点坐标。

基于MTLAB的3节点三角形单元求解平面弹性问题

本文部分内容来自网络整理，本司不为其真实性负责，如有异议或侵权请及时联系，本司将立即删除！== 本文为word格式，下载后可方便编辑和修改！ ==深知兵真爱兵,打牢官兵思想篇一：对战士的爱更深沉些对战士的爱更深沉些、对战士关心更长远些一直以来，总队医院党委始终把促进战士成长、成才作为提高部队战斗力的重要举措。

针对战士整体素质高但个体差异较大、所学知识多但弄懂搞精的少的实际，医院党委想方设法为战士学习成才创造条件，尽力满足战士求知成才愿望。

使医院开展的“深知兵、真爱兵、带好兵”活动迸发出更大活力。

一、在“深知”上下功夫。

知兵是爱兵的基础。

总队医院党委不仅要求机关干部和带兵干部在了解大量的基本数据，掌握思想变化的“活情况”的同时，而且深层次了解战士能干什么，给每名战士“量体裁衣”，为每个战士建立成才档案，每年在新战士下连后，总队医院都会开展“走近士兵、研究士兵”活动。

主动摸清战士需求，要对战士的长远发展负责，引导他们把个人追求与部队建设有机统一起来，在岗位上成才，在成才中进步。

为战士合理制定成才计划，使战士的合理愿望切合实际。

设身处地的关爱，拉近了官兵心与心之间的距离，战士有话愿向干部讲，有事愿靠组织办，有苦愿向骨干诉，有乐愿与大家享的氛围已经形成，上下和谐，官兵关系顺畅。

二、在“真爱”上做文章。

总队医院积极为战士成才创造条件，医院党委投资100多万元为中队建立电脑网络学习室，与驻地4家大型企业建立了战士成才协作机制。

总队医院充分运用医院自身人才技术队伍力量，对战士进行计算机操作、医学理疗、超声波技术、放射技术等专项技能培训。

为了使培训取得实效，使战士的计算机技能得到社会认可，医院还主动邀请、积极协调地方人力资源部门专业人员到部队授课，并对参加培训的战士进行了职业技能鉴定考核。

目前，参加各类培训的战士陆续如愿获得各种职业资格证书。

同时，政治处鼓励战士参加军地自考、函授学习，并借助各种优势，邀请地方院校教授来部队作专题知识讲座，为战士学习成才搭建平台，调动了战士学习成才积极性。

matlab有限元三角形单元编程

matlab有限元三角形单元编程
在MATLAB中进行有限元分析，可以使用其自带的有限元分析工具箱（FEATool）进行编程。

以下是一个简单的例子，演示如何使用三角形单元进行有限元分析：
1. 打开MATLAB，进入FEATool环境。

2. 创建新的有限元模型。

选择“Model”菜单下的“New Model”选项，进入“Model Builder”界面。

3. 在“Model Builder”界面中，选择“2D Triangle”单元类型，并在绘图区域中绘制出三角形网格。

4. 在“Model Builder”界面中，设置材料属性、边界条件和载荷等参数。

5. 运行有限元分析。

选择“Model”菜单下的“Solve”选项，进行有限元求解。

6. 查看结果。

选择“Model”菜单下的“Results”选项，可以查看位移、应力、应变等结果。

以上是一个简单的例子，演示了如何使用三角形单元进行有限元分析。

在实际应用中，还需要根据具体问题进行详细的建模和计算。

三角形单元数值积分 matlab

三角形单元数值积分 matlab
在Matlab中进行三角形单元数值积分可以通过使用内置的函数
来实现。

一种常用的方法是使用`integral`函数来进行数值积分。

假设我们有一个三角形单元的函数f(x)，我们可以使用以下步骤来
进行数值积分：
步骤1，定义三角形单元的函数f(x)。

这可能涉及到使用三角
形的顶点坐标和函数值来定义一个插值函数。

步骤2：使用`integral`函数对定义的函数f(x)进行数值积分。

例如，如果我们的函数是f(x)，我们可以使用以下命令来进行数值
积分：
matlab.
integral(@(x) f(x), a, b)。

其中a和b是积分的下限和上限。

步骤3，根据需要，可以使用不同的数值积分方法，例如
'auto'（自动选择方法）、'tiled'（瓦片方法）或者
'ArrayValued'（对数组进行积分）等。

另外，如果需要对三角形单元进行数值积分，也可以考虑使用`trapz`函数进行梯形数值积分。

这可以通过将三角形边界上的点作为离散数据点来实现。

需要注意的是，在使用Matlab进行三角形单元数值积分时，需要确保对积分区域进行合适的离散化，以便进行数值计算。

同时，也需要考虑数值积分的精度和误差控制，以确保得到准确的积分结果。

总之，Matlab提供了丰富的数值积分函数和方法，可以方便地对三角形单元进行数值积分，用户可以根据具体情况选择合适的方法来进行数值积分计算。

6节点三角形单元matlab

6节点三角形单元matlab在MATLAB中，可以使用有限元分析工具箱来创建和操作三角形有限元网格以及进行相应的计算。

下面我将从几个角度来介绍在MATLAB中如何处理6节点三角形单元。

1. 创建6节点三角形单元网格：在MATLAB中，可以使用 PDE 模型创建器 App 或者命令行函数来创建6节点三角形单元网格。

首先，你可以使用命令行函数如下创建一个简单的6节点三角形单元网格：matlab.model = createpde();geometryFromEdges(model,@()gunit);generateMesh(model,'Hmax',0.1);pdeplot(model);这段代码创建了一个 PDE 模型，使用默认的单位几何体，然后生成了一个包含6节点三角形单元的网格，并进行了绘制。

2. 定义材料属性和边界条件：一旦创建了网格，你可以定义材料属性和边界条件。

例如，你可以使用以下代码定义一个材料属性和施加边界条件：matlab.specifyCoefficients(model,'m',0,'d',1,'c',1,'a',0,'f',1);applyBoundaryCondition(model,'dirichlet','Edge',1:model.Geo metry.NumEdges,'u',0);3. 求解和后处理：接下来，可以使用 solve 函数求解模型，并使用结果进行后处理。

例如，可以使用以下代码求解模型并绘制结果：matlab.result = solvepde(model);pdeplot(model,'XYData',result.NodalSolution);这将绘制出6节点三角形单元的有限元解。

总结：在MATLAB中处理6节点三角形单元，首先需要创建网格，然后定义材料属性和边界条件，最后求解模型并进行后处理。

matlab有限元三角形单元程序

matlab有限元三角形单元程序以下是一个简单的 MATLAB 有限元三角形单元程序的示例：```matlab% 定义模型参数E = 1000; % 弹性模量nu = 0.3; % 泊松比thickness = 1; % 板的厚度% 定义节点坐标nodes = [0, 0; 1, 0; 1, 1; 0, 1; 0.5, 0.5];% 定义单元连接关系connectivity = [1, 2, 5; 2, 3, 5; 3, 4, 5; 4, 1, 5];% 计算总节点数和总单元数numNodes = size(nodes, 1);numElements = size(connectivity, 1);% 初始化全局刚度矩阵和载荷向量K = zeros(numNodes);F = zeros(numNodes, 1);% 循环遍历每个单元for i = 1:numElements% 查找当前单元的节点编号nodesIndex = connectivity(i, :);% 根据节点编号从全局坐标矩阵中取出节点坐标coordinates = nodes(nodesIndex, :);% 计算当前单元的局部刚度矩阵localK = calculateLocalStiffness(E, nu, thickness, coordinates);% 组装局部刚度矩阵到全局刚度矩阵中K(nodesIndex, nodesIndex) = K(nodesIndex, nodesIndex) + localK;% 计算当前单元的局部载荷向量localF = calculateLocalLoad(thickness, coordinates);% 组装局部载荷向量到全局载荷向量中F(nodesIndex) = F(nodesIndex) + localF;end% 边界条件：节点1固定K(1, :) = 0;K(1, 1) = 1;F(1) = 0;% 解线性方程组U = K \ F;% 输出位移结果disp('节点位移：');disp(U);% 计算应力结果stress = calculateStress(E, nu, thickness, nodes, connectivity, U);% 输出应力结果disp('节点应力：');disp(stress);% 计算局部刚度矩阵的函数function localK = calculateLocalStiffness(E, nu, thickness, coordinates)% 计算单元的雅可比矩阵J = (1/2) * [coordinates(2,1)-coordinates(1,1), coordinates(3,1)-coordinates(1,1);coordinates(2,2)-coordinates(1,2), coordinates(3,2)-coordinates(1,2)];% 计算雅可比矩阵的逆矩阵invJ = inv(J);% 计算单元刚度矩阵B = invJ * [-1, 1, 0; -1, 0, 1];D = (E/(1-nu^2)) * [1, nu, 0; nu, 1, 0; 0, 0, (1-nu)/2]; localK = thickness * abs(det(J)) * (B' * D * B);end% 计算局部载荷向量的函数function localF = calculateLocalLoad(thickness, coordinates) localF = zeros(3, 1);midPoint = [sum(coordinates(:,1))/3,sum(coordinates(:,2))/3];localF(3) = thickness * 1 * det([coordinates(1,:); coordinates(2,:); midPoint]);end% 计算各节点应力的函数function stress = calculateStress(E, nu, thickness, nodes, connectivity, U)stress = zeros(size(nodes, 1), 3);for i = 1:size(connectivity, 1)nodesIndex = connectivity(i, :);coordinates = nodes(nodesIndex, :);Ke = calculateLocalStiffness(E, nu, thickness, coordinates);Ue = U(nodesIndex);stress(nodesIndex, :) = stress(nodesIndex, :) + (Ke * Ue)';endstress = stress / thickness;end```这个程序实现了一个简单的平面三角形单元有限元分析，包括定义节点坐标和单元连接关系、计算全局刚度矩阵和载荷向量、施加边界条件、解线性方程组、计算节点位移和应力等。

3结点三角形单元有限元程序MATLAB语言

3结点三角形单元有限元程序（MATLAB语言）该程序包括以下6个部分：1．主程序tri_fem：用于数据的录入和其他程序的调用；2．总刚程序Kf：计算结构的总体刚度；3．各结点位移求解程序xf：求解各结点的位移；4．线性方程组求解程序Jordan：Gauss-Jordan法求解非约束结点的位移；5．应力应变程序ss：由各结点位移求解各单元内的三个结点的应力stress和应变strain；6．数据录入程序input：录入材料、几何尺寸、单元编号和结点编号、位移约束和已知载荷等。

以课本P25页例2.2为例，其input程序为function [E,v,t,EN,ecode,NN,node,RN,RC,PN,PC]=input()E=2.1e11; v=1/3; t=1; %杨氏模量Pa，泊松比，厚度EN=2; %单元数ecode=[3 1 2; %单元编号单元1 3-1-2;单元2 1-3-41 3 4];NN=4; %结点数node=[0 0; %各结点坐标2 0;2 1;0 1];RN=2; %被约束的位移数RC=[1 4]; %被约束的结点PN=2; %有载荷的结点数%PC（1）表示有载荷的结点，PC(2)表示各结点的力，PC（3）表示载荷方向，0为x方向，1为y方向PC=[2 3;-1/2 -1/2;1 1]; %结点2、3分别有y负方向上-1/2N的力作用在matlab环境下，输入则程序运行结果如下：该程序求解的结点位移结果和结点应力结果与课本给出的结果一致。

附录：%%-------------平面三角形单元有限元法---------------------%%function [x strain stress]=tri_fem()[E,v,t,EN,ecode,NN,node,RN,RC,PN,PC]=input; %调入已定材料、几何尺寸以及单元和结点编号及约束和载荷分布[n m]=size(ecode);if EN~=nerror('Wrong elementnumber EN or wrong elementcode ecode!');return;end[n m]=size(node);if NN~=nerror('Wrong nodenumber NN or wrong node-coordinate node!');return;ende=zeros(EN,6);A=zeros(EN,1); %面积for i=1:ENe(i,:)=[node(ecode(i,1),:),node(ecode(i,2),:),node(ecode(i,3),:)]; %各三角形单元的节点坐标D=[1,node(ecode(i,1),:)1,node(ecode(i,2),:)1,node(ecode(i,3),:)];A(i,1)=1/2*det(D);end%% 形成单元参数b=zeros(EN,3);c=zeros(EN,3); %各单元参数初始化for i=1:ENb(i,1)=e(i,4)-e(i,6); b(i,2)=e(i,6)-e(i,2); b(i,3)=e(i,2)-e(i,4);c(i,1)=-e(i,3)+e(i,5); c(i,2)=-e(i,5)+e(i,1); c(i,3)=-e(i,1)+e(i,3);end%% 求得总刚，并引入约束和载荷求得各结点位移K=Kf(E,v,t,EN,ecode,NN,A,b,c); %调用函数Kf，求得结构的总体刚度矩阵x=xf(NN,RN,RC,PN,PC,K); %调用函数xf,求得各结点位移%% 求解应力应变[strain stress]=ss(E,v,EN,ecode,A,b,c,x);%% 单元刚度矩阵与结构刚度矩阵function K=Kf(E,v,t,EN,ecode,NN,A,b,c)Ke=zeros(6,6); %单元的刚度矩阵，初始为6*6阶零矩阵K=zeros(NN*2,NN*2); %结构的总体刚度矩阵，初始为零矩阵for m=1:1:EN %m为单元号for i=1:1:3for j=1:1:3Ke(2*i-1,2*j-1)=b(m,i)*b(m,j)+(1-v)*c(m,i)*c(m,j)/2;Ke(2*i-1,2*j)=v*c(m,i)*b(m,j)+(1-v)*b(m,i)*c(m,j)/2;Ke(2*i,2*j-1)=v*b(m,i)*c(m,j)+(1-v)*c(m,i)*b(m,j)/2;Ke(2*i,2*j)=c(m,i)*c(m,j)+(1-v)*b(m,i)*b(m,j)/2;endendKe=E*t/(4*(1-v^2)*A(EN)).*Ke; %获得单元m的刚度矩阵Kb=mat2cell(Ke,ones(1,3)*2,ones(1,3)*2); %将单元矩阵Ke分为3*3块set1=ones(1,NN)*2;Ka=mat2cell(zeros(NN*2,NN*2),set1,set1); %将总刚K分为NN*NN块for i=1:1:3for j=1:1:3Ka(ecode(m,i),ecode(m,j))=Kb(i,j); %各单元刚度矩阵整体编号,并叠加endendK=K+cell2mat(Ka);end%分块矩阵K合成一个矩阵K%% 引入位移向量和右端项function x=xf(NN,RN,RC,PN,PC,K)x=ones(NN*2,1); %位移初始为0向量for i=1:RNx(RC(i)*2-1)=0;x(RC(i)*2)=0;end %被约束结点位移为0%%----------------引入已知结点载荷-------------%px=zeros(NN*2,1); %载荷初始为0向量for i=1:PNif PC(3,i)==1px(PC(1,i)*2)=PC(2,i);else if PC(3,i)==0px(PC(1,i)*2-1)=PC(2,i);endendend%%----------------引入已知结点载荷-------------%%%----------------求解非约束结点的位移X-------------%set1=ones(1,NN)*2;Ka=mat2cell(K,set1,set1);pxa=mat2cell(px,set1,1);AN=zeros(2*(NN-RN),2*(NN-RN));ANa=mat2cell(AN,ones(1,NN-RN)*2,ones(1,NN-RN)*2);bn=zeros(2*(NN-RN),1);bna=mat2cell(bn,ones(1,NN-RN)*2,1);BN=zeros(2*RN,2*(NN-RN));BNa=mat2cell(BN,ones(1,RN)*2,ones(1,NN-RN)*2);m=1;for i=1:1:NNif x(2*i)==1M(m)=i;m=m+1;endendfor i=1:1:NN-RNfor j=1:1:NN-RNANa(i,j)=Ka(M(i),M(j));bna(i,1)=pxa(M(i),1);endendfor i=1:RNfor j=1:NN-RNBNa(i,j)=Ka(RC(i),M(j));endendAN=cell2mat(ANa);bn=cell2mat(bna);BN=cell2mat(BNa);X=Jordan(AN,bn); %利用Gauss-Jordan法求解非约束结点的位移X %%----------------求解非约束结点的位移X-------------%%----------------由X可得所有结点位移x-------------%BN=BN*X;m=1; n=1;for i=1:1:NNif x(2*i)==1x(2*i-1)=X(m);x(2*i)=X(m+1);m=m+2;else if x(2*i)==0px(2*i-1)=BN(n);px(2*i)=BN(n+1);n=n+2;endendend%% 列主元Jordan消去法将系数矩阵化成对角矩阵求解方程组的数值解function x=Jordan(A,b)%开始计算，赋初值[n,m]=size(A);x=zeros(n,1);for k=1:n%选主元max1=0;for i=k:nif abs(A(i,k))>max1max1=abs(A(i,k));r=i;endend%交换两行if r>kfor j=k:nz=A(k,j); A(k,j)=A(r,j); A(r,j)=z;endz=b(k); b(k)=b(r); b(r)=z;end%消元计算b(k)=b(k)/A(k,k);for j=k+1:nA(k,j)=A(k,j)/A(k,k);endfor i=1:nif i~=kfor j=k+1:nA(i,j)=A(i,j)-A(i,k)*A(k,j);endb(i)=b(i)-A(i,k)*b(k);endendend%输出xfor i=1:nx(i)=b(i);end%% 求解应力应变function [strain stress]=ss(E,v,EN,ecode,A,b,c,x)strain=zeros(3,EN); %应变初始为0矩阵stress=zeros(3,EN); %应力初始为0矩阵D=E/(1-v^2)*[1 v 0;v 1 0;0 0 (1-v)/2];for m=1:1:ENB=[b(m,1) 0 b(m,2) 0 b(m,3) 0;0 c(m,1) 0 c(m,2) 0 c(m,3);c(m,1) b(m,1) c(m,2) b(m,2) c(m,3) b(m,3)];B=B/2/A(m,1); %应变矩阵S=D*B; %应力矩阵X=[x(2*ecode(m,1)-1),x(2*ecode(m,1)),x(2*ecode(m,2)-1),x(2*ecode(m,2)),x(2*ecode(m,3)-1),x (2*ecode(m,3))]';strain(:,m)=B*X;stress(:,m)=S*X;end。

清华大学弹性力学-大作业matlab代码

f1=int(N',m,-1,1); f=1/(4*a*b)*2*int(f1,n,-1,1); f=subs(f); %重积分过程求单元结点荷载 f %初始给定 K 为零阵
K=zeros((a+1)*(b+1),(a+1)*(b+1)); for x1=1:(1+a) for y1=1:(1+b) for x2=1:(1+a) for y2=1:(1+b) s1=(y1-1)*(a+1)+x1; s2=(y2-1)*(a+1)+x2; if x1==x2 & y1==y2 if x1==1
T=round(T); K=K(T,T); F=F(T); U=K\F; U(1,:) %A 点位移
end if abs(x1-x2)==1 & abs(y1-y2)==1 if (x1-x2)*(y1-y2)==1 K(s1,s2)=k(2,4); else K(s1,s2)=k(1,3); end end end end end end K; %以上循环生成 K 矩阵 %如果需要可以通过此处查看 K
1/3]; A\b 3.a×b 个矩形单元 %a×b 个矩形单元 clc,clear all a=input('单元格长数'); a=round(a); b=input('单元格宽数'); b=round(b); x=0:1/a:1; y=0:1/b:1; A=[1 -1 1 -1 1111 1 1 -1 -1 1 -1 -1 1]; P=inv(A); syms m n; %m,n 相当于ξ和η N=[1,m,n,m*n]*P; B=[2*a*(P(2,:)+n*P(4,:)) 2*b*(P(3,:)+m*P(4,:))]; D=eye(2); k1=int(B'*D*B,m,-1,1); k=1/(4*a*b)*int(k1,n,-1,1); k=subs(k); %重积分过程求单元刚度矩阵 k

三角形有限单元法程序

global gNdt gElt gNTu
hold on dFau=1.0; %变形图放大系数 %绘制变形前网格 trisurf(gElt,gNdt(:,1),gNdt(:,2),zeros(size(gNdt,1),1)) view(2); axis equal; axis off; axis tight; %获取变形后数据 DDisp=gNdt+gNTu*dFau; pause(3.0); trisurf(gElt,DDisp(:,1),DDisp(:,2),zeros(size(DDisp,1),1)); view(2); axis equal; axis off; axis tight;
－练习题1 图(a)所示的分析物体，根据对称型取1/4作为分析模型，如图 (b)所示。划分为4个三角形单元，6个结点。 y
y
2N/m
2
该问题可以简化为平面应力问题求解，假设弹性模量为E=1，泊松比为0 。该问题的解析解答为：
i
m
j
%计算刚度矩阵
Ke=Bm'*De*Bm*A;
－单元刚度矩阵的组集
1 2 3 4 5 6 7 8 9 10 11 12
2 3 4 5 6 7 8 9 10 11 12 1
2
i
* * * * * * * * * * * * * * * *
gMAT-----材料链表
gElt-------各单元的结点编号
gEInfo----各单元的材料号 gNdt------各结点的坐标 gBct-------边界条件(位移和力)
gKA-------组集后的总体刚度矩阵
gFA--------组集后的总体荷载列阵 gNTu------求解后得到的结点位移

单元刚度矩阵(等参元)MATLAB编程

《有限元法》实验报告专业班级力学（实验）1601 姓名田诗豪学号 1603020210 提交日期 2019.4.24实验一（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=0.5*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,0.5,1.0];%****输入材料参数矩阵(弹性模量，泊松比，壁厚)****[B3,S3,K3]=ele_mat_tri3(xy3,mat)三、算例分析●算例1：如图1所示三角形单元，结点坐标为1(0,0)，2(5,2)，3(1,4)，弹性模量为200GPa，泊松比为0.35、厚度为0.5m。

matlab程序画几何图形并进行网格划分

编写程序，对下图几何区域按有限元法进行网格划分，确定单元中结点数目及位置，并对单元和结点分别进行编号。

二、程序思路基于MATLAB软件平台，首先画出几何图形，然后对几何区域按照有限元法思想进行三角形网格剖分，最后按顺序依次输出单元及结点编号和坐标。

三、编写源程序：clear;%第一步：画出几何区域%为按逆时针顺序绘制几何区域x_up=[3.5,0];y_up=[2,2];%上固体边界数组x_l=[0,0];y_l=[2,0];%入流左边界数组x_low=[0,2.5];y_low=[0,0];%下固体边界数组n=9;%圆弧等分变量nc=linspace(pi,pi/2,n);r=1;%圆弧半径变量rx_c=r*cos(c)+3.5;y_c=r*sin(c);%三角函数表示圆柱左上部分x_r=[3.5,3.5];y_r=[1,2];%出流右边界数组x=[x_low,x_c,x_r,x_up,x_l];%整体边界X坐标数组y=[y_low,y_c,y_r,y_up,y_l];%整体边界Y坐标数组plot(x,y,'m');%绘制几何区域xlabel('x轴');ylabel('y轴');title('几何区域');hold onl=k*(pi/2)*r;x_up(3)=3.5-l;y_up(3)=2;cut_point=[x_up(3),y_up(3)];x_up(3)=x_up(2);y_up(3)=y_up(2);x_up(2)=cut_point(1);y_up(2)=cut_point(2);%将上边界点坐标按顺序重排x_cutline=[x_up(2),x_low(2)];y_cutline=[y_up(2),y_low(2)];%分割线数组plot(x_cutline,y_cutline);%绘制分割线n1=6;%上边界左半部分n1等分变量x_upcp1=linspace(x_up(2),x_up(3),n1+1);%上边界左半部分n1等分y_upcp1=linspace(y_up(2),y_up(3),n1+1);x_lowcp1=linspace(x_low(2),x_low(1),n1+1);%下边界n1等分y_lowcp1=linspace(y_low(2),y_low(1),n1+1);x_upcp2=linspace(x_up(1),x_up(2),n+1);%上边界右半部分n等分y_upcp2=linspace(y_up(1),y_up(2),n+1);x_upper=[x_upcp2,x_upcp1];%上固体边界坐标数组y_upper=[y_upcp2,y_upcp1];c1=linspace(pi*3/2,pi,n+1);x_c1=cos(c1)+3.5;y_c1=abs(sin(c1));x_lower=[x_c1,x_lowcp1];y_lower=[y_c1,y_lowcp1];for i=1:n1+n+2%画竖线plot([x_upper(i)x_lower(i)],[y_upper(i)y_lower(i)],'m');endm=7;for i=n+n1+2:-1:1plot([x_upper(i)x_lower(i)],[y_upper(i)y_lower(i)],'m');%绘制竖线x1(i,:)=linspace(x_upper(i),x_lower(i),m+1);%竖线m等分y1(i,:)=linspace(y_upper(i),y_lower(i),m+1);endfor i=1:mfor j=n+n1+2:-1:2plot([x1(j,i)x1(j-1,i)],[y1(j,i)y1(j-1,i)],'m');%plot([x1(j,i)x1(j-1,i+1)],[y1(j,i)y1(j-1,i+1)],'m');endendfor j=1:n+n1for i=1:mp1=i+(j-1)*(m+1);p2=i+1+j*(m+1);p3=i+j*(m+1);t=2*i+2*(j-1)*m;p4=i+(j-1)*(m+1);p5=i+1+(j-1)*(m+1);p6=i+1+j*(m+1);e(t-1,:)=[p1,p2,p3];%整体编号与局部编号关系e(t,:)=[p4,p5,p6];endend四、程序各步运行结果：具体程序如上所写第一步：画出几何区域。

算例分析(1)

解：（1）结构的离散化与编号对该结构进行离散，单元编号及节点编号如图4-20(b)所示，即有二个3节点三角形单元。

载荷F 按静力等效原则向节点1、节点2移置等效。

matlab编译平面有限元计算

matlab编译平面有限元计算编译平面有限元计算是一种常用的数值计算方法，可以用于求解各种复杂的工程问题。

在本文中，我们将介绍如何使用MATLAB编写平面有限元计算程序，并通过一个实例来说明其应用。

让我们了解一下有限元方法的基本原理。

有限元方法是一种将连续体划分为有限个单元，根据物理方程和边界条件，在每个单元上建立离散方程组，最终求解得到整个连续体的近似解的方法。

在平面问题中，连续体被划分为三角形或四边形单元。

MATLAB是一种功能强大的数值计算软件，它提供了丰富的工具箱和函数，可以方便地进行有限元计算。

下面我们将以求解平面弹性力学问题为例，介绍如何使用MATLAB编写平面有限元计算程序。

我们需要定义问题的几何信息、边界条件和材料参数。

在MATLAB中，可以通过定义节点坐标、单元连接关系、边界条件和材料参数来描述问题。

节点坐标可以用一个矩阵表示，其中每一行代表一个节点的坐标。

单元连接关系可以用一个矩阵表示，其中每一行代表一个单元的节点编号。

边界条件可以用一个向量表示，其中每个元素代表一个节点的边界条件。

材料参数可以用一个矩阵表示，其中每一行代表一个单元的材料参数。

接下来，我们需要建立有限元离散方程组。

在平面弹性力学问题中，离散方程组可以通过组装单元刚度矩阵和外力向量得到。

单元刚度矩阵可以通过单元的几何信息和材料参数计算得到。

外力向量可以根据边界条件和载荷信息计算得到。

最终，我们可以通过求解离散方程组得到节点的位移解。

我们可以根据节点的位移解计算出单元的应力和应变。

在平面弹性力学问题中，应力可以通过单元的几何信息和位移解计算得到。

应变可以通过应力和材料参数计算得到。

通过计算应力和应变，我们可以评估结构的稳定性和性能。

MATLAB编译平面有限元计算程序是一种强大的工具，可以用于求解各种复杂的工程问题。

通过定义几何信息、边界条件和材料参数，建立有限元离散方程组，求解得到节点的位移解，最终计算出单元的应力和应变。

三节点三角形单元matlab编程

1,'%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,NPOIN])'
% 结点坐标数组
FORCE=fscanf(FP1,'%f',[3,NFORCE])' % 结点力数组
FIXED=fscanf(FP1,'%d',[3,inf])'
% 约束数组
ASTIF=ASSEMBLE(NPOIN,NELEM,YOUNG,POISS,THICK,COORD,LNODS) % 生成总刚
%总刚度矩阵
function ASTIF=ASSEMBLE(NPOIN,NELEM,YOUNG,POISS,THICK,COORD,LNODS)
% 张成特定大小总刚矩阵，并赋初值 0
ASTIF(1:2*NPOIN,1:2*NPOIN)=0;
% 组装总刚度矩阵
for i=1:NELEM
% 计算单刚
[ESTIF,SMATX]=FORMESTIF(i,COORD,LNODS,YOUNG,POISS,THICK)
总结点数, 单元数,
受约束点数, 结点力数 ,
第一组：
NPOIN=fscanf(FP1,'%d',1),NELEM=fscanf(FP1,'%d',1),
NFORCE=fscanf(FP1,'%d',1),NVFIX=fscanf(FP1,'%d',1)

matlab三角形单元扩展编程

matlab三角形单元扩展编程在MATLAB 中，扩展三角形单元的编程通常涉及到对现有的三角形单元进行修改或添加新的功能。

下面是一个简单的例子，演示如何创建一个用于扩展三角形单元的MATLAB 类。

假设你想要创建一个扩展的三角形单元类，其中包含了额外的属性和方法。

以下是一个可能的实现示例：```matlabclassdef ExtendedTriangle < handlepropertiesNodes % 顶点坐标Area % 三角形面积endmethods% 构造函数function obj = ExtendedTriangle(nodes)if nargin > 0obj.Nodes = nodes;obj.calculateArea();endend% 计算三角形面积function calculateArea(obj)% 在这里实现计算三角形面积的逻辑% 这可能涉及到对顶点坐标的操作% 这里仅作为示例，假设nodes 是一个3x2 的矩阵obj.Area = polyarea(obj.Nodes(:,1), obj.Nodes(:,2));end% 自定义的其他方法可以在这里添加endend```在这个示例中，`ExtendedTriangle` 类具有`Nodes` 属性表示三角形的顶点坐标，并且有一个`Area` 属性表示三角形的面积。

构造函数负责初始化这些属性，而`calculateArea` 方法用于计算三角形的面积。

你可以根据实际需求添加更多的属性和方法。

使用这个类的示例：```matlab% 创建一个三角形对象triangleNodes = [0, 0; 1, 0; 0, 1];triangle = ExtendedTriangle(triangleNodes);% 访问属性disp('Triangle Nodes:');disp(triangle.Nodes);disp(['Triangle Area: ', num2str(triangle.Area)]);```这只是一个简单的例子，你可以根据实际需要扩展和修改这个类。

平面梁单元MATLAB有限元程序

平面梁单元MATLAB有限元程序function []=beam2n()clear allclose allclearclose%------------------------ BEAM2 ---------------------------disp('========================================'); disp(' PROGRAM BEAM2 '); disp(' Beam Bending Analysis '); disp(' T.R.Chandrupatla and A.D.Belegundu '); disp('========================================');InputData;Bandwidth;Stiffness;ModifyForBC;BandSolver;ReactionCalc;Output;%------------------------ function InputData ---------------------------function [] = InputData();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTglobal TITLE FILE1 FILE2global LINP LOUTglobal NQdisp(blanks(1));FILE1 = input('Input Data File Name ','s'); LINP = fopen(FILE1,'r');FILE2 = input('Output Data File Name ','s'); LOUT = fopen(FILE2,'w');DUMMY = fgets(LINP);TITLE = fgets(LINP);DUMMY = fgets(LINP);TMP = str2num(fgets(LINP));[NN, NE, NM, NDIM, NEN, NDN] =deal(TMP(1),TMP(2),TMP(3),TMP(4),TMP(5),TMP(6));NQ = NDN * NN;DUMMY = fgets(LINP);TMP = str2num(fgets(LINP));[ND, NL, NMPC]= deal(TMP(1),TMP(2),TMP(3));NPR=1; % E%----- Coordinates -----DUMMY = fgets(LINP);for I=1:NNTMP = str2num(fgets(LINP));[N, X(N,:)]=deal(TMP(1),TMP(2:1+NDIM)); end%----- Connectivity -----DUMMY = fgets(LINP);for I=1:NETMP = str2num(fgets(LINP));[N,NOC(N,:), MAT(N), SMI(N)] = ...deal(TMP(1),TMP(2:1+NEN), TMP(2+NEN), TMP(3+NEN));end%----- Specified Displacements ----- DUMMY = fgets(LINP); for I=1:NDTMP = str2num(fgets(LINP));[NU(I,:),U(I,:)] = deal(TMP(1), TMP(2)); end%----- Component Loads -----DUMMY = fgets(LINP);F = zeros(NQ,1);for I=1:NLTMP = str2num(fgets(LINP));[N,F(N)]=deal(TMP(1),TMP(2)); end%----- Material Properties ----- DUMMY = fgets(LINP);for I=1:NMTMP = str2num(fgets(LINP));[N, PM(N,:)] = deal(TMP(1), TMP(2:NPR+1));end%----- Multi-point Constraints B1*Qi+B2*Qj=B0if NMPC > 0DUMMY = fgets(LINP);for I=1:NMPCTMP = str2num(fgets(LINP));[BT(I,1), MPC(I,1), BT(I,2), MPC(I,2), BT(I,3)] = ...deal(TMP(1),TMP(2),TMP(3),TMP(4),TMP(5));endendfclose(LINP);%------------------------ function Bandwidth ---------------------------function []=Bandwidth();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT STRESS REACT global CNSTglobal TITLE FILE1 FILE2global LINP LOUT%----- Bandwidth Evaluation ----- NBW = 0;for N=1:NENABS = NDN*(abs(NOC(N, 1) - NOC(N, 2)) + 1);if (NBW < NABS)NBW = NABS;endendfor I=1:NMPCNABS = abs(MPC(I, 1) - MPC(I, 2)) + 1;if (NBW < NABS)NBW = NABS;endenddisp(sprintf('Bandwidth = %d', NBW));%------------------------ function Stiffness ---------------------------function []=Stiffness();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBW global X NOC F AREA MAT SMI S global PM NU U MPC BT REACTglobal CNSTglobal NQ%----- Global Stiffness Matrix S = zeros(NQ,NBW);for N = 1:NEdisp(sprintf('Forming Stiffness Matrix of Element %d', N));%-------- Element Stiffness and Temperature Load -----N1 = NOC(N, 1);N2 = NOC(N, 2);M = MAT(N);EL = abs(X(N1) - X(N2));EIL = PM(M, 1) * SMI(N) / EL^3;SE(1, 1) = 12 * EIL;SE(1, 2) = EIL * 6 * EL;SE(1, 3) = -12 * EIL;SE(1, 4) = EIL * 6 * EL;SE(2, 1) = SE(1, 2);SE(2, 2) = EIL * 4 * EL * EL;SE(2, 3) = -EIL * 6 * EL;SE(2, 4) = EIL * 2 * EL * EL;SE(3, 1) = SE(1, 3);SE(3, 2) = SE(2, 3);SE(3, 3) = EIL * 12;SE(3, 4) = -EIL * 6 * EL;SE(4, 1) = SE(1, 4);SE(4, 2) = SE(2, 4);SE(4, 3) = SE(3, 4);SE(4, 4) = EIL * 4 * EL * EL;disp('.... Placing in Global Locations'); for II = 1:NENNRT = NDN * (NOC(N, II) - 1);for IT = 1:NDNNR = NRT + IT;I = NDN * (II - 1) + IT;for JJ = 1:NENNCT = NDN * (NOC(N, JJ) - 1);for JT = 1:NDNJ = NDN * (JJ - 1) + JT;NC = NCT + JT - NR + 1;if (NC > 0)S(NR, NC) = S(NR, NC) + SE(I, J);endendendendendend%------------------------ function ModifyForBC ---------------------------function []=ModifyForBC();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTglobal NQ%----- Decide Penalty Parameter CNST ----- CNST = 0;for I = 1:NQif CNST < S(I, 1); CNST = S(I, 1); end endCNST = CNST * 10000;%----- Modify for Boundary Conditions ----- % --- Displacement BC ---for I = 1:NDN = NU(I);S(N, 1) = S(N, 1) + CNST;F(N) = F(N) + CNST * U(I);end%--- Multi-point Constraints ---for I = 1:NMPCI1 = MPC(I, 1);I2 = MPC(I, 2);S(I1, 1) = S(I1, 1) + CNST * BT(I, 1) * BT(I, 1);S(I2, 1) = S(I2, 1) + CNST * BT(I, 2) * BT(I, 2);IR = I1;if IR > I2; IR = I2; endIC = abs(I2 - I1) + 1;S(IR, IC) = S(IR, IC) + CNST * BT(I, 1) * BT(I, 2);F(I1) = F(I1) + CNST * BT(I, 1) * BT(I, 3);F(I2) = F(I2) + CNST * BT(I, 2) * BT(I, 3); end%------------------------ function BandSolver ---------------------------function []=BandSolver();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTglobal NQ%----- Equation Solving using Band Solver ----- disp('Solving using Band Solver(bansol.m)');[F] = bansol(NQ,NBW,S,F);%------------------------ function ReactionCalc ---------------------------function []=ReactionCalc();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTfor I = 1:NDN = NU(I);REACT(I) = CNST * (U(I) - F(N));end%------------------------ function Output ---------------------------function []=Output();global NN NE NM NDIM NEN NDNglobal ND NL NCH NPR NMPC NBWglobal X NOC F AREA MAT SMI Sglobal PM NU U MPC BT REACTglobal CNSTglobal TITLE FILE1 FILE2global LINP LOUTdisp(sprintf('Output for Input Data from file %s\n',FILE1)); fprintf(LOUT,'Output for Input Data from file %s\n',FILE1);disp(TITLE);fprintf(LOUT,'%s\n',TITLE);disp(' Node# X-Displ Rotation');fprintf(LOUT,' Node# X-Displ Rotation\n'); I=[1:NN]';% print a matrixdisp(sprintf(' %4d %15.4E %15.4E\n',[I,F(2*I-1),F(2*I)]')); fprintf(LOUT,' %4d %15.4E %15.4E\n',[I,F(2*I-1),F(2*I)]');%----- Reaction Calculation -----disp(sprintf(' DOF# Reaction'));fprintf(LOUT,' DOF# Reaction\n');for I = 1:NDN = NU(I);R = CNST * (U(I) - F(N));disp(sprintf(' %4d %15.4E',N,REACT(I)));fprintf(LOUT,' %4d %15.4E\n',N,REACT(I)); endfclose(LOUT);disp(sprintf('The Results are available in the text file %s', FILE2));。

单元刚度矩阵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，泊松比为、厚度为。

平面三角形单元常应变单元matlab程序的编制

平面三角形单元常应变单元matlaｂ程序的编制
———————————————————————————————— 作者：
———————————————————————————————— 日期：
ﻩ
三角形常应变单元程序的编制与使用
有限元法是求解微分方程边值问题的一种通用数值方法,该方法是一种基于变分法(或变分里兹法）而发展起来的求解微分方程的数值计算方法，以计算机为手段，采用分片近似，进而逼近整体的研究思想求解物理问题。
有限元法中三节点三角形分析结构的步骤如下:
1）整理原始数据,如材料性质、荷载条件、约束条件等,离散结构并进行单元编码、结点编码、结点位移编码、选取坐标系。
2）单元分析,建立单元刚度矩阵。
3）整体分析，建立总刚矩阵。
4）建立整体结构的等效。
POISＳ=fsｃanf（FＰ1，'%f'，1）%泊松比
TＨICＫ=ｆscａnf（ＦＰ1,'%d',1)％厚度
LＮODS=fscanf（ＦＰ1,'%ｄ',［3,NＥLEM］)'%单元定义数组(单元结点号）
功能:运用有限元法中三角形常应变单元解平面问题的计算主程序。
基本思想：单元结点按右手法则顺序编号。
荷载类型：可计算结点荷载。
说明:主程序的作用是通过赋值语句、读取和写入文件、函数调用等完成算法的全过程，即实现程序流程图的程序表达。
%---－--－-----－---－-----－--－---－－---－--－-----－-－-－－－－-－-－－-－-----------------－----－----－----－----------
glｏｂａｌFＯRCＥFIXＥＤ
globaｌＢMATＸＤMＡTＸSＭATＸAＲＥA