有限元基礎編程（終結篇）——C3D8單元程序編制

易公子

2022年7月24日 15:48

瀏覽：7310 評論：42 收藏：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]