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

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

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

因此，本文首先重點(diǎn)講解的溫度作用產(chǎn)生的節(jié)點(diǎn)等效溫度荷載有限元列式的推導(dǎo)，六面體單元?jiǎng)偠染仃嚒⒀鸥鞅染仃嚒?yīng)變矩陣有限元列式的推導(dǎo)，溫度應(yīng)力應(yīng)變的后處理有限元列式，此外還包括上述有限元列式對應(yīng)的matlab代碼實(shí)現(xiàn)過程。

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

假定通過傳熱分析我們可以得到，相對于原來狀態(tài)溫度升高了，對于各向同性材料，溫度升高會(huì)產(chǎn)生一個(gè)均勻的應(yīng)變（通俗地講，就是熱脹冷縮原理，物體會(huì)發(fā)生膨脹），應(yīng)變大小與材料的線膨脹系數(shù)有關(guān)，表示材料在單位溫度升高所引起的長度的變化值，一般情況下，在一定范圍內(nèi)可以假定線膨脹系數(shù)為常數(shù)。若物體能夠自由變形，則由溫度變化引起的應(yīng)變不會(huì)產(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)系代入到有限元分析列式中，利用虛功原理進(jìn)行推導(dǎo)。其中虛應(yīng)力和虛應(yīng)變場函數(shù)的有限元列式如下式

? ? ? ? ? ?（3）

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

? ? ? （4）

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

? ??（5）

因?yàn)樘撐灰拼嬖谌我庑裕虼丝梢猿艄剑?）中的虛位移，進(jìn)而得到

? ? ?（6）

其中Ke為剛度矩陣，表達(dá)式為

? ? ? ?（7）

為單元節(jié)點(diǎn)載荷，表達(dá)式為

? ? ??（8）

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

??? ? ? （9）

因此，對于溫度應(yīng)力問題，核心就是求解等效溫度載荷。因?yàn)楣街邪薆矩陣，即應(yīng)變矩陣，對于六面體單元應(yīng)變矩陣的推導(dǎo)過程如下圖所示。

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

基于圖1所示的應(yīng)變矩陣的推導(dǎo)過程，公式（9）中等效溫度載荷對應(yīng)的matlab代碼如下所示，此外本段代碼同樣包含了單元?jiǎng)偠染仃嚨木幊虒?shí)現(xiàn)，因?yàn)榈刃囟容d荷和單元?jiǎng)偠染仃嚨那蠼夤街芯蠦矩陣和D矩陣參與。

for X=1:2
for Y=1:2
for Z=1:2
GP1=GaussCoordinate(X); GP2=GaussCoordinate(Y); GP3=GaussCoordinate(Z); %高斯點(diǎn)坐標(biāo)
% 計(jì)算形函數(shù)對總體坐標(biāo)的導(dǎo)數(shù)（NDerivative）及雅可比矩陣行列式（JacobiDET）
[~,NDerivative, JacobiDET] = ShapeFunction([GP1 GP2 GP3], ElementNodeCoordinate);
Coefficient=GaussWeight(X)*GaussWeight(Y)*GaussWeight(Z)*JacobiDET;
%計(jì)算B矩陣  利用形函數(shù)對總體坐標(biāo)的導(dǎo)數(shù)（NDerivative）對B進(jìn)行計(jì)算
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ù)導(dǎo)數(shù)的ShapeFunction函數(shù)的matlab代碼如下所示：

function [N,NDerivative,JacobiDET] = ShapeFunction(GaussPoint,ElementNodeCoordinate)
%等參元坐標(biāo)  每一列代表一個(gè)點(diǎn)的坐標(biāo)
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ù)對局部坐標(biāo)的導(dǎo)數(shù)矩陣3*8
%計(jì)算形函數(shù)及形函數(shù)對局部坐標(biāo)導(dǎo)數(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;%計(jì)算雅可比矩陣
JacobiDET = det(Jacobi);%計(jì)算雅可比行列式
JacobiINV=inv(Jacobi);%對雅可比行列式求逆
NDerivative=JacobiINV*ParentNDerivative;%利用雅可比行列式的逆計(jì)算形函數(shù)對結(jié)構(gòu)坐標(biāo)的導(dǎo)數(shù)3*8
end

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

InterpolationMatrix=zeros(8,8);%求解節(jié)點(diǎn)應(yīng)力應(yīng)變的插值矩陣
%循環(huán)高斯點(diǎ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;
% 計(jì)算局部坐標(biāo)下的形函數(shù)及形函數(shù)導(dǎo)數(shù)
[N,NDerivative, ~] = ShapeFunction([E1 E2 E3], ElementNodeCoordinate);
ElementNodeDisplacement=U(ElementNodeDOF);
ElementNodeDisplacement=reshape(ElementNodeDisplacement,Dof,8);
% 計(jì)算高斯點(diǎn)應(yīng)變 GausspointStrain3_3 3*3的應(yīng)變矩陣
GausspointStrain3_3=ElementNodeDisplacement*NDerivative’;%3*8  8*3
%把高斯點(diǎn)應(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)]';
% 計(jì)算高斯點(diǎn)應(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é)點(diǎn)應(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)力問題的整個(gè)求解過程，以一個(gè)兩端固定約束的四棱柱為例，其在不同單位熱應(yīng)變的作用下的求解結(jié)果如下圖所示

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