有限元基础编程（终结篇）——C3D8单元程序编制
有限元基础编程（终结篇）——C3D8单元程序编制
本篇推文是有限元基础编程的终结篇，讲述C3D8单元的程序编制及实现。主要内容有：C3D8单元理论基础、便于编程的“乘大数法”处理边界条件、编制程序注意事项、云图绘制函数、INP文件读取函数、Abaqus仿真对比等，内容量大，慢慢食用~
特别声明：程序框架采用了吉林大学左文杰老师的脚本文件，计算单元刚度的核心计算程序仍延续我们以往编制程序的风格。
理论基础
与Q4单元理论基础相同，唯一的区别就是：每个节点的自由度由2变成了3，代码具体变化看Ke函数和C3D8_cal_B函数的变化，理论部分可参考有限元基础编程——Q4单元。
边界条件处理（乘大数法）
理论：
对于边界条件
考察其第
由于
即
优点：
1. 既可以处理
的情形，又可以处理 的情形； 2. 原整体刚度矩阵求解规模不变，不需要重新排序；
3. 保持整体刚度矩阵的对称性。
相关代码：
% 乘大数法施加位移约束
BigNumber=1e8;
ConstraintsNumber=size(Constraints,1);
FixedDof=Dof*(Constraints(:,1)-1)+Constraints(:,2); %被约束的自由度编号(列向量)
for i=1:ConstraintsNumber
K(FixedDof(i),FixedDof(i))=K(FixedDof(i),FixedDof(i))*BigNumber;
Force(FixedDof(i))=Constraints(i,3)*K(FixedDof(i),FixedDof(i));
end程序编制
模型描述：一悬臂梁，左端固定，三个方向自由度均受约束，右端下侧受到均布荷载作用。
图2 悬臂梁示意图
主程序代码
%****************************************************
% ***：易木木响叮当
%****************************************************
function Main()
%读取inp文件获得节点坐标信息Nodes及单元信息Elements
[Nodes, Elements] = Readmesh( 'BC_Force.inp' );
Forces=[601 2 -100;607 2 -100;613 2 -100;619 2 -100;625 2 -100];
ConNumber=1:30;
Constraints=zeros(size(ConNumber,2)*3,3);
for i=1:size(ConNumber,2)
Constraints(3*i-2:3*i,:)=[ConNumber(i) 1 0;ConNumber(i) 2 0;ConNumber(i) 3 0;];
end
E=210000; %弹性模量
u=0.3; %泊松比
U =StaticsSolver(E,u,Forces,Constraints,Nodes,Elements);
% 输出结果
OutputTXT = fopen('Results.txt','w'); %打开一个可写文件，用于写入计算结果
OutputResults(OutputTXT,Nodes,Elements,U)%调用输出结果文件
fclose(OutputTXT);
edit('Results.txt')
end以上程序采用了吉林大学左文杰老师编制的脚本文件，需要注意以下几个方面：
1. 程序中使用
Readmesh函数（仅需调用即可）读取INP文件中节点、单元信息；2. 节点外载荷和边界条件信息需要在INP文件中寻找，较为麻烦，方便起见可采用ANSYS建模，导出节点、单元、载荷、边界条件信息，相应的APDL命令流如下：
nwrite,node,txt,
ewrite,element,txt,
Dlist,ALL
Flist,ALL1. 云图绘制采用
PlotContour函数，本程序仅仅绘制出网格图、位移云图，因为研究对象是线弹性体，位移是第一求解未知量，最为准确，应力是通过应力矩阵变换而得，线弹性状态很简单，不是求解核心，所以本程序仅展示位移云图，有关应力云图显示详情请观看一阶六面体(C3D8)的线弹性有限元编程[1]
Abaqus对比
云图对比
图3 Abaqus-U2 位移云图
图4 Matlab-U2 位移云图
Matlab最大U2值为-0.0312，Abaqus最大U2值为-0.0320，结果相差2.5%
引用链接
[1] 一阶六面体(C3D8)的线弹性有限元编程: https://www.bilibili.com/video/BV1KU4y1U7Mz?spm_id_from=333.999.0.0&vd_source=21e509b445f316000076bb38c1f19a05
