有限元基础编程(终结篇)——C3D8单元程序编制

有限元基础编程(终结篇)——C3D8单元程序编制

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

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

理论基础

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

边界条件处理(乘大数法)

理论

对于边界条件     的情形,可将刚度方程进行如下操作:

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

考察其第     行的等价性: 

   

由于    (   是一个很大的数),则上式变为

   

即    

优点

  1. 1. 既可以处理     的情形,又可以处理     的情形;

  2. 2. 原整体刚度矩阵求解规模不变,不需要重新排序;

  3. 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(i1 0;ConNumber(i2 0;ConNumber(i3 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. 1. 程序中使用Readmesh函数(仅需调用即可)读取INP文件中节点、单元信息;

  2. 2. 节点外载荷和边界条件信息需要在INP文件中寻找,较为麻烦,方便起见可采用ANSYS建模,导出节点、单元、载荷、边界条件信息,相应的APDL命令流如下:

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

Abaqus对比

云图对比

有限元基础编程(终结篇)——C3D8单元程序编制的图3

图3 Abaqus-U2 位移云图

有限元基础编程(终结篇)——C3D8单元程序编制的图4

图4 Matlab-U2 位移云

Matlab最大U2值为-0.0312,Abaqus最大U2值为-0.0320,结果相差2.5%

代码文件:在公众号:易木木响叮当,后台回复C3D8,即可自动获取~

引用链接

[1] 一阶六面体(C3D8)的线弹性有限元编程: https://www.bilibili.com/video/BV1KU4y1U7Mz?spm_id_from=333.999.0.0&vd_source=21e509b445f316000076bb38c1f19a05

(35条)
默认 最新
编程难度很大
评论 1 点赞 2
回复
Dei
评论 点赞
厉害的易公子,给力,谢谢分享,支持一波
评论 1 点赞 1
回复
感谢
评论 点赞

查看更多评论 >

点赞 65 评论 92 收藏 37
关注