有限元基础编程（终结篇）——C3D8单元程序编制

易公子 2022年7月24日浏览：3984 评论：92 收藏：37

有限元基础编程（终结篇）——C3D8单元程序编制

本篇推文是有限元基础编程的终结篇，讲述C3D8单元的程序编制及实现。主要内容有：C3D8单元理论基础、便于编程的“乘大数法”处理边界条件、编制程序注意事项、云图绘制函数、INP文件读取函数、Abaqus仿真对比等，内容量大，慢慢食用~

特别声明：程序框架采用了吉林大学左文杰老师的脚本文件，计算单元刚度的核心计算程序仍延续我们以往编制程序的风格。代码文件获取方式详见文末。

理论基础

与Q4单元理论基础相同，唯一的区别就是：每个节点的自由度由2变成了3，代码具体变化看Ke函数和C3D8_cal_B函数的变化，理论部分可参考有限元基础编程——Q4单元。

边界条件处理（乘大数法）

理论：

对于边界条件 $\bar{q}r=\bar{u}$ 的情形，可将刚度方程进行如下操作：

有限元基础编程（终结篇）——C3D8单元程序编制的图1

考察其第 $r$ 行的等价性：

$k{r1}q_1+k_{r2}q_2+\cdots +\alpha k_{rr}\bar{q}r+\cdots +k{rn}q_n=\alpha k_{rr}\bar{u}$

由于 $\alpha k_{rr}\gg k_n$ ( $\alpha$ 是一个很大的数)，则上式变为

$\alpha k_{rr}\bar{q}r\approx \alpha k{rr}\bar{u}$

即 $\bar{q}_r\approx \bar{u}$

优点：

1. 既可以处理 $\bar{q}_r=0$ 的情形，又可以处理 $\bar{q}_r=\bar{u}$ 的情形；
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

程序编制

模型描述：一悬臂梁，左端固定，三个方向自由度均受约束，右端下侧受到均布荷载作用。

有限元基础编程（终结篇）——C3D8单元程序编制的图2

图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,ALL

1. 云图绘制采用PlotContour函数，本程序仅仅绘制出网格图、位移云图，因为研究对象是线弹性体，位移是第一求解未知量，最为准确，应力是通过应力矩阵变换而得，线弹性状态很简单，不是求解核心，所以本程序仅展示位移云图，有关应力云图显示详情请观看一阶六面体(C3D8)的线弹性有限元编程^[1]