关于平面应力问题有限元法——三角形常应变单元Matlab求解器

根据有限元中三角形常应变单元的理论知识，利用Matlab工具，编制了本求解器，用以求解平面应力问题，并将所求结果与Ansys软件所求结果及弹性理学中的理论值进行了比较，效果较好。

% “三角形单元“求解“平面应力问题“——处理器部分
% 说明：在运行本程序之前，请做好以下准备工作-预处理
% 1、预处理
% 1.0 实常数定义
%     弹性模量 EX
%     泊松比   PRXY
%     板厚     t
% 注意：以上各量单位封闭，自成体系    
% 1.1 三角形单元位移模式 u=a1+a2*x+a3*y;
%                       v=a4+a5*x+a6*y;
% 1.2 弹性体剖分，划分三角形单元   1、2...T （T 表示三角形的个数）
%                对节点进行标号   1、2...n （n 表示节点的个数）
%                给出节点坐标信息 Nodes[i1]（i=1..n） 表示第i个节点的“x“坐标
%                                Nodes[i2]（i=1..n） 表示第i个节点的“y“坐标
%                写出单元连接信息 Element（ki） （k=1..T） 表示第“k“个三角形的“i“节点
%                                Element（kj） （k=1..T） 表示第“k“个三角形的“j“节点
%                                Element（km） （k=1..T） 表示第“k“个三角形的“m“节点
% 1.3 外载荷等效节点载荷 Load
%                       Load（ixy）（i=1..n） 表示i节点的等效载荷沿 “x“、“y“方向的分力

% 2、处理器
% 2.1 计算单元刚度矩阵   
%     （1）、[K]=（66i）（i=1..T）表示第i个三角形的刚度矩阵 
%     （2）、[K]=t*AREA（i）*[B]‘[D][B]
%          其中 [D]为平面应力问题中的物理矩阵
%               [B]为几何矩阵，可通过下式求出
% 2.2 组装形成整体刚度矩阵  [k]=（2*n2*n）即大小为2倍的节点数的方阵
% 2.3 引入位移边界条件  
%     Hold（2*节点个数1）=1 表示该位移被限制， Hold（2*节点个数11）=0 表示为自由位移
% 2.4 模型的求解
%     （1）消元法    即“划行划列法“
%     （2）置大数法
%     说明，本程序中采用的 前一种方法

% 3、后处理
% 3.1 求应力         Thigma     最大应力 Max_Thigma
% 3.2 绘制变形后网络
% 3.3 绘制出应力场图
% 3.4 导出图表、网页等

%界面初始化
clc;
clear;

%加载初始信息
load Indata.mat;
% ~~EX_弹性模量
% ~~Element_单元连接信息
% ~~Hold_节点位移约束
% ~~Nodes_节点坐标
% ~~PRXY_泊松比
% ~~t_板厚

% 计算物理矩阵 D（33）;
D=EX/（1-PRXY^2）*[1 PRXY 0;PRXY 1 0;0 0 （1-PRXY）/2];

% 计算每个三角形的面积 AREA
  AREA=zeros（size（Element1）1）;      % size（Element2） 测量划分三角形的个数
  for i=1: size（Element1）
      AREA（i）=1/2*det（...
            [ 1Nodes（Element（i1）1）Nodes（Element（i1）2）; ...
              1Nodes（Element（i2）1）Nodes（Element（i2）2）; ...
              1Nodes（Element（i3）1）Nodes（Element（i3）2）]）; 
  end    
% 定义单元刚度矩阵 大小为6*6
  K=zeros（66size（Element1））;
% 定义单元几何矩阵
  B=zeros（36size（Element1））;
  for i=1:size（Element1）
      delta=...
      [ 1Nodes（Element（i1）1）Nodes（Element（i1）2）; ...
        1Nodes（Element（i2）1）Nodes（Element（i2）2）; ...
        1Nodes（Element（i3）1）Nodes（Element（i3）2）]; 
      for j=1:3
          %其中 Math（deltaji）表示求delta 矩阵j行i列的代数余子式
          B（12*j-1i）=Math（deltaj2）*1/2/AREA（i）;
          B（32*ji）=Math（deltaj2）*1/2/AREA（i）;
          B（22*ji）=Math（deltaj3）*1/2/AREA（i）;
          B（32*j-1i）=Math（deltaj3）*1/2/AREA（i）;
      end
  end
% 计算单元刚度矩阵
  for i=1:size（Element1）
      K（::i）=t*AREA（i）*B（::i）‘*D*B（::i）;
  end
% 叠加得整体刚度矩阵
Node=0;                    % 标记节点序号
% 定义整体刚度矩阵
k=zeros（2*size（Nodes1）2*size（Nodes1））;
%单元刚度矩阵的投放
for i=1:size（Element1）
    for j=1:3
    Node=Element（ij）;
    k（2*Node-1:2*Node2*Node-1:2*Node）=k（2*Node-1:2*Node2*Node-1:2*Node）+K（2*j-1:2*j2*j-1:2*ji）;
    end
    k（2*Element（i1）-1:2*Element（i1）2*Element（i2）-1:2*Element（i2））=...
        k（2*Element（i1）-1:2*Element（i1）2*Element（i2）-1:2*Element（i2））+K（1:23:4i）;
    k（2*Element（i1）-1:2*Element（i1

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件       8275  2010-01-09 19:43  In.txt

     文件       1290  2010-01-10 15:43  Indata.mat

     文件       5487  2010-01-10 16:15  Main.m

     文件       5488  2010-01-10 15:54  Main.asv

     文件        289  2009-12-26 12:35  Math.m

     文件        230  2009-12-26 11:31  Math.asv

     文件       4217  2010-01-09 19:41  Out.txt

     文件     142177  2010-01-09 19:28  实验报告.pdf

----------- ---------  ---------- -----  ----

               167453                    8