基于六面體單元熱應(yīng)力問題的Matlab有限元編程求解

導讀：上一篇《彈性地基梁matlab有限元編程，以雙排樁支護結(jié)構(gòu)計算為例》引起了Matlab有限元編程學習者的共鳴。今天我想和大家討論一下熱應(yīng)力問題（六面體單元）matlab有限元編程問題。

本案例主要介紹熱應(yīng)力問題的matlab有限元編程及原理。熱應(yīng)力也叫溫度應(yīng)力，后文提到的熱應(yīng)變也叫溫度應(yīng)變。這里需要與傳熱問題的有限元分析進行區(qū)分：可以認為傳熱分析是進行熱應(yīng)力分析的前提條件，通過傳熱分析來確定溫度場，在獲得溫度場的基礎(chǔ)上，計算所產(chǎn)生的熱應(yīng)力。所以本文為闡述熱應(yīng)力問題前提是假定溫度場是已知的，我們并不關(guān)心溫度場是如何得到的，那是熱傳導要解決的問題，在這個已知的溫度場作用下，由于熱脹冷縮會產(chǎn)生響應(yīng)的熱應(yīng)變進而產(chǎn)生熱應(yīng)力，如何求解這個熱應(yīng)力，這就是我們這次課程要解決的問題。

如圖1所示，本案例的分析對象為一個兩端固定約束的四棱柱受不同方向的熱應(yīng)變作用下產(chǎn)生的應(yīng)力應(yīng)變，在用有限元方法求解熱應(yīng)力問題時，溫度作用可以等效為一個節(jié)點載荷進行求解，因此熱應(yīng)力問題本質(zhì)上還是一個固體力學問題的有限元求解，而我們剛才提到的傳熱問題的有限元求解其實是在求解傅里葉傳熱偏微分方程，其與固體力學偏微分方程是完全不同的兩套理論。

因此，本文首先重點講解的溫度作用產(chǎn)生的節(jié)點等效溫度荷載有限元列式的推導，六面體單元剛度矩陣、雅各比矩陣、應(yīng)變矩陣有限元列式的推導，溫度應(yīng)力應(yīng)變的后處理有限元列式，此外還包括上述有限元列式對應(yīng)的matlab代碼實現(xiàn)過程。

圖1 兩端固定約束的四棱柱受單位熱應(yīng)變作用下的應(yīng)力

假定通過傳熱分析我們可以得到，相對于原來狀態(tài)溫度升高了，對于各向同性材料，溫度升高會產(chǎn)生一個均勻的應(yīng)變（通俗地講，就是熱脹冷縮原理，物體會發(fā)生膨脹），應(yīng)變大小與材料的線膨脹系數(shù)有關(guān)，表示材料在單位溫度升高所引起的長度的變化值，一般情況下，在一定范圍內(nèi)可以假定線膨脹系數(shù)為常數(shù)。若物體能夠自由變形，則由溫度變化引起的應(yīng)變不會產(chǎn)生應(yīng)力。溫度應(yīng)變一般被表示為初應(yīng)變的形式

? ??（1）

對應(yīng)的應(yīng)力-應(yīng)變關(guān)系為

? ? ? ? （2）

接下來將公式（2）的應(yīng)力應(yīng)變關(guān)系代入到有限元分析列式中，利用虛功原理進行推導。其中虛應(yīng)力和虛應(yīng)變場函數(shù)的有限元列式如下式

? ? ? ? ? ?（3）

將公式（2）利用虛功原理，可得到下述虛功方程

? ? ? （4）

將公式（3）代入公式（4）所示的虛功方程中，并對其進行簡化處理，可以得到

? ??（5）

因為虛位移存在任意性，因此可以除掉公式（5）中的虛位移，進而得到

? ? ?（6）

其中Ke為剛度矩陣，表達式為

? ? ? ?（7）

為單元節(jié)點載荷，表達式為

? ? ??（8）

為等效溫度載荷，是由熱應(yīng)力引起的節(jié)點處等效溫度載荷，表達式為

??? ? ? （9）

因此，對于溫度應(yīng)力問題，核心就是求解等效溫度載荷。因為公式中包含了B矩陣，即應(yīng)變矩陣，對于六面體單元應(yīng)變矩陣的推導過程如下圖所示。

圖1 應(yīng)變矩陣與剛度矩陣的推導公式

基于圖1所示的應(yīng)變矩陣的推導過程，公式（9）中等效溫度載荷對應(yīng)的matlab代碼如下所示，此外本段代碼同樣包含了單元剛度矩陣的編程實現(xiàn)，因為等效溫度載荷和單元剛度矩陣的求解公式中均有B矩陣和D矩陣參與。

for X=1:2
for Y=1:2
for Z=1:2
GP1=GaussCoordinate(X); GP2=GaussCoordinate(Y); GP3=GaussCoordinate(Z); %高斯點坐標
% 計算形函數(shù)對總體坐標的導數(shù)（NDerivative）及雅可比矩陣行列式（JacobiDET）
[~,NDerivative, JacobiDET] = ShapeFunction([GP1 GP2 GP3], ElementNodeCoordinate);
Coefficient=GaussWeight(X)*GaussWeight(Y)*GaussWeight(Z)*JacobiDET;
%計算B矩陣  利用形函數(shù)對總體坐標的導數(shù)（NDerivative）對B進行計算
B=zeros(6,24);
for I=1:8
COL=(I-1)*3+1:(I-1)*3+3;
B(:,COL)=[NDerivative(1,I) 0         0;
0         NDerivative(2,I) 0;
0         0         NDerivative(3,I);
NDerivative(2,I) NDerivative(1,I) 0;
0         NDerivative(3,I) NDerivative(2,I);
NDerivative(3,I) 0         NDerivative(1,I)];
end
Ke=Ke+Coefficient*B'*D*B;  %疊加剛度陣
P_dT=P_dT+Coefficient*B'*D*epsilon0; %等效溫度荷載
end
end
end

上述代碼中，求解雅各比矩陣行列式和形函數(shù)導數(shù)的ShapeFunction函數(shù)的matlab代碼如下所示：

function [N,NDerivative,JacobiDET] = ShapeFunction(GaussPoint,ElementNodeCoordinate)
%等參元坐標  每一列代表一個點的坐標
ParentNodes=[-1 1  1 -1 -1  1  1 -1;
-1 -1  1  1 -1 -1  1  1;
-1 -1 -1 -1  1  1  1  1];
N=zeros(8,1); %初始化形函數(shù)矩陣8*1
ParentNDerivative=zeros(3,8);%初始化形函數(shù)對局部坐標的導數(shù)矩陣3*8
%計算形函數(shù)及形函數(shù)對局部坐標導數(shù)
for I=1:8
XPoint = ParentNodes(1,I);
YPoint = ParentNodes(2,I);
ZPoint = ParentNodes(3,I);
ShapePart = [1+GaussPoint(1)*XPoint 1+GaussPoint(2)*YPoint 1+GaussPoint(3)*ZPoint];  %定義中間變量
N(I) = 0.125*ShapePart(1)*ShapePart(2)*ShapePart(3);
ParentNDerivative(1,I) = 0.125*XPoint*ShapePart(2)*ShapePart(3);
ParentNDerivative(2,I) = 0.125*YPoint*ShapePart(1)*ShapePart(3);
ParentNDerivative(3,I) = 0.125*ZPoint*ShapePart(1)*ShapePart(2);
end
Jacobi = ParentNDerivative*ElementNodeCoordinate;%計算雅可比矩陣
JacobiDET = det(Jacobi);%計算雅可比行列式
JacobiINV=inv(Jacobi);%對雅可比行列式求逆
NDerivative=JacobiINV*ParentNDerivative;%利用雅可比行列式的逆計算形函數(shù)對結(jié)構(gòu)坐標的導數(shù)3*8
end

單元剛度矩陣和等效溫度載荷求得之后，將其裝配到整體剛度矩陣和荷載向量中通過高斯消去法實現(xiàn)節(jié)點位移的求解，求解之后的節(jié)點位移按照下圖中的后處理流程和對應(yīng)的公式，即可求解得到高斯積分點處的應(yīng)力應(yīng)變值，需要注意的是如公式（2）所示單元應(yīng)力的求解公式中包含了溫度應(yīng)變；下一步再通過插值函數(shù)矩陣將高斯積分點處的應(yīng)力應(yīng)變插值到每個節(jié)點再進行磨平處理，即可得到節(jié)點的應(yīng)力應(yīng)變。具體matlab實現(xiàn)代碼如下所示：

InterpolationMatrix=zeros(8,8);%求解節(jié)點應(yīng)力應(yīng)變的插值矩陣
%循環(huán)高斯點
for X=1:2
for Y=1:2
for Z=1:2  
E1=GaussCoordinate(X); E2=GaussCoordinate(Y); E3=GaussCoordinate(Z);
GaussPointNumber = GaussPointNumber + 1;
% 計算局部坐標下的形函數(shù)及形函數(shù)導數(shù)
[N,NDerivative, ~] = ShapeFunction([E1 E2 E3], ElementNodeCoordinate);
ElementNodeDisplacement=U(ElementNodeDOF);
ElementNodeDisplacement=reshape(ElementNodeDisplacement,Dof,8);
% 計算高斯點應(yīng)變 GausspointStrain3_3 3*3的應(yīng)變矩陣
GausspointStrain3_3=ElementNodeDisplacement*NDerivative’;%3*8  8*3
%把高斯點應(yīng)變矩陣改寫成6*1
GausspointStrain=[GausspointStrain3_3(1,1) GausspointStrain3_3(2,2) GausspointStrain3_3(3,3)…  
GausspointStrain3_3(1,2)+GausspointStrain3_3(2,1)…
GausspointStrain3_3(2,3)+GausspointStrain3_3(3,2) …
GausspointStrain3_3(1,3)+GausspointStrain3_3(3,1)]';
% 計算高斯點應(yīng)力
GausspointStress = D*GausspointStrain-D*epsilon0;%總應(yīng)變-溫度應(yīng)變
GaussStrain(:,GaussPointNumber)=GausspointStrain;
GaussStress(:,GaussPointNumber)=GausspointStress;
InterpolationMatrix(K,:)=N;
K=K+1;
end
end
end
%求得節(jié)點應(yīng)力應(yīng)變
Temp1=InterpolationMatrix(GaussStrain(1:6,GaussPointNumber-7:GaussPointNumber)');
NodeStrain(1:6,INODE+1:INODE+8)=Temp1';
Temp2=InterpolationMatrix(GaussStress(1:6,GaussPointNumber-7:GaussPointNumber)');
NodeStress(1:6,INODE+1:INODE+8)=Temp2';

上述便是熱應(yīng)力問題的整個求解過程，以一個兩端固定約束的四棱柱為例，其在不同單位熱應(yīng)變的作用下的求解結(jié)果如下圖所示

以上就是筆者本期要分享的內(nèi)容。