平面三角形单元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中，我们可以使用矩阵变换的方法来实现三角形的旋转。
本文将介绍如何使用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。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
%变量说明
%NPOIN NELEM NVFIX NFORCE NNODE
%总结点数，单元数，受约束边界点数，结点力数, 单元结点数
%COORD LNODS YOUNG POISS THICK %结构结点坐标数组，单元定义数组，弹性模量，泊松比，厚度
%FORCE FIXED BMA TX DMATX SMATX %节点力数组，约束信息数组，单元应变矩阵，单元弹性矩阵，单元应力矩阵%AREA HK ASLOD ASDISP FP1
%单元面积，总体刚度矩阵，总体荷载向量，结点位移向量，数据文件指针format short e
clear
FP1=fopen('C:\Users\Administrator\Desktop\input.txt','rt');
NPION=fscanf(FP1,'%d',1); %结点个数（结点编码总数）NELEM=fscanf(FP1,'%d',1); %单元个数（单元编码总数）NFORCE=fscanf(FP1,'%d',1); %结点荷载个数
NVFIX=fscanf(FP1,'%d',1); %受约束边界点数
YOUNG=fscanf(FP1,'%e',1); %弹性模量
POISS=fscanf(FP1,'%f',1); %泊松比
THICK=fscanf(FP1,'%d',1); %厚度
LNODS=fscanf(FP1,'%d',[3,NELEM])'; %单元定义数组（单元结点号）COORD=fscanf(FP1,'%f',[2,NPION])'; %结点坐标数组
FORCE=fscanf(FP1,'%f',[3,NFORCE])'; %结点力数组
FIXED=fscanf(FP1,'%d',[3,inf])'; %约束信息数组
%引用所需的全局变量
%global NPION NELEM COORD LNODS YOUNG POISS
%global BMA TX DMATX SMA TX AREA
%生成弹性矩阵D
a=YOUNG/(1-POISS^2);
DMATX(1,1)=1*a;
DMATX(1,2)=POISS*a;
DMATX(2,1)=POISS*a;
DMATX(2,2)=1*a;
DMATX(3,3)=(1-POISS)*a/2;
for i=1:NELEM; %i为当前所计算的单元号
%计算当前单元的面积
AREA=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;
end
%生成应变矩阵B
for j=0:2
b(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);
end
BMA TX=[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*AREA);
SMA TX=DMATX*BMATX; %求应力矩阵S=D*B
% 所引用的全局变量：
%global NPION NELEM NVFIX LNODS ASTIF THICK
%global BMA TX SMA TX AREA FIXED
HK=seros(2*NPION,2*NPION); %张成特定大小总刚矩阵并置0
for i=1:NELEM
EK=BMA TX'*SMATX*THICK*AREA; %求解单元刚度矩阵
a=LNODS(i,:); %临时向量,用来记录当前单元的节点编号
for j=1:3
for k=1:3
HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)=HK((a(j)*2-1):a(j)*2,(a(k)*2-1):a(k)*2)+EK(j*2-1:j*2,k* 2-1:k*2);
%跟据节点编号对应关系将单元刚度分块叠加到总刚矩阵中
end
end
end
%将约束信息加入总刚（置0置1法）
NUM=1; %计数器(当前已分析的节点数)
i=0; %计数器(当前已处理的约束数)
tmp(NVFIX)=0; %临时存被约束的序号
while i<NVFIX
for j=-1:0
if FIXED(NUM,j+3)==1 %若发现约束
i=i+1; %计数器自增
tmp(i)=FIXED(NUM)*2+j; %求约束序号
end
end
NUM=NUM+1;
end
for i=1:NVFIX
HK(1:2*NPION,tmp(i))=0; %将一约束序号处的总刚列向量清0
HK(tmp(i),1:2*NPION)=0; %.将一约束序号处的总刚行向量清0
HK(tmp(i),tmp(i))=1; %将行列交叉处的元素置为1
end
%所需引用的全局变量
%global ASLOD NPION NFORCE FORCE
ASLOD(1:2*NPION)=0; %张成特定大小的0向量
for i=1:NFORCE
ASLOD((FORCE(i,1)*2-1):FORCE(i,1)*2)=FORCE(i,2:3);
end
%将有约束处的荷载置为0
NUM=1;
i=0;
tmp(NVFIX)=0;
while i<NVFIX
for j=-1:0
if FIXED(NUM,j+3)==1
i=i+1;
tmp(i)=FIXED(NUM)*2+j;
end
end
NUM=NUM+1;
end
for i=1:NVFIX
ASLOD(tmp(i))=0;
end
ASDISP=HK\ASLOD'; %计算结点位移向量
%所引用的全局变量
%global NELEM NPION SMA TX ASDISP LNODS
ELEDISP(1:6)=0; %当前单元的结点位移向量for i=1:NELEM
for j=1:3
ELEDISP(j*2-1:j*2)=ASDISP(LNODS(i,j)*2-1:LNODS(i,j)*2);
end
STRESS=SMA TX*ELEDISP'; %求内力end。