点击文尾阅读原文试看
首发 | 仿真秀App
本案例主要介绍热应力问题的matlab有限元编程及原理。热应力也叫温度应力，后文提到的热应变也叫温度应变。这里需要与传热问题的有限元分析进行区分：可以认为传热分析是进行热应力分析的前提条件，通过传热分析来确定温度场，在获得温度场的基础上，计算所产生的热应力。所以本文为阐述热应力问题前提是假定温度场是已知的，我们并不关心温度场是如何得到的，那是热传导要解决的问题，在这个已知的温度场作用下，由于热胀冷缩会产生响应的热应变进而产生热应力，如何求解这个热应力，这就是我们这次课程要解决的问题。
如图1所示，本案例的分析对象为一个两端固定约束的四棱柱受不同方向的热应变作用下产生的应力应变，在用有限元方法求解热应力问题时，温度作用可以等效为一个节点载荷进行求解，因此热应力问题本质上还是一个固体力学问题的有限元求解，而我们刚才提到的传热问题的有限元求解其实是在求解傅里叶传热偏微分方程，其与固体力学偏微分方程是完全不同的两套理论。
因此，本文首先重点讲解的温度作用产生的节点等效温度荷载有限元列式的推导，六面体单元刚度矩阵、雅各比矩阵、应变矩阵有限元列式的推导，温度应力应变的后处理有限元列式，此外还包括上述有限元列式对应的matlab代码实现过程。
图1 两端固定约束的四棱柱受单位热应变作用下的应力
，对于各向同性材料，温度升高
会产生一个均匀的应变（通俗地讲，就是热胀冷缩原理，物体会发生膨胀），应变大小与材料的线膨胀系数
有关，
表示材料在单位温度升高所引起的长度的变化值，一般情况下，在一定范围内可以假定线膨胀系数为常数。若物体能够自由变形，则由温度变化引起的应变不会产生应力。温度应变一般被表示为初应变的形式
（1）
对应的应力-应变关系为
（2）
接下来将公式（2）的应力应变关系代入到有限元分析列式中，利用虚功原理进行推导。其中虚应力和虚应变场函数的有限元列式如下式
（3）
将公式（2）利用虚功原理，可得到下述虚功方程
（4）
将公式（3）代入公式（4）所示的虚功方程中，并对其进行简化处理，可以得到
（5）
因为虚位移存在任意性，因此可以除掉公式（5）中的虚位移，进而得到
（6）
其中Ke为刚度矩阵，表达式为
（7）
为单元节点载荷，表达式为
（8）
为等效温度载荷，是由热应力引起的节点处等效温度载荷，表达式为
（9）
因此，对于温度应力问题，核心就是求解等效温度载荷。因为公式中包含了B矩阵，即应变矩阵，对于六面体单元应变矩阵的推导过程如下图所示。
图1 应变矩阵与刚度矩阵的推导公式
基于图1所示的应变矩阵的推导过程，公式（9）中等效温度载荷对应的matlab代码如下所示，此外本段代码同样包含了单元刚度矩阵的编程实现，因为等效温度载荷和单元刚度矩阵的求解公式中均有B矩阵和D矩阵参与。
for X=1:2for Y=1:2for Z=1:2GP1=GaussCoordinate(X); GP2=GaussCoordinate(Y); GP3=GaussCoordinate(Z); %高斯点坐标% 计算形函数对总体坐标的导数（NDerivative）及雅可比矩阵行列式（JacobiDET）[~,NDerivative, JacobiDET] = ShapeFunction([GP1 GP2 GP3], ElementNodeCoordinate);Coefficient=GaussWeight(X)*GaussWeight(Y)*GaussWeight(Z)*JacobiDET;%计算B矩阵 利用形函数对总体坐标的导数（NDerivative）对B进行计算B=zeros(6,24);for I=1:8COL=(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)];endKe=Ke+Coefficient*B'*D*B; %叠加刚度阵P_dT=P_dT+Coefficient*B'*D*epsilon0; %等效温度荷载endendend
上述代码中，求解雅各比矩阵行列式和形函数导数的ShapeFunction函数的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); %初始化形函数矩阵8*1ParentNDerivative=zeros(3,8);%初始化形函数对局部坐标的导数矩阵3*8%计算形函数及形函数对局部坐标导数for I=1:8XPoint = 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);endJacobi = ParentNDerivative*ElementNodeCoordinate;%计算雅可比矩阵JacobiDET = det(Jacobi);%计算雅可比行列式JacobiINV=inv(Jacobi);%对雅可比行列式求逆NDerivative=JacobiINV*ParentNDerivative;%利用雅可比行列式的逆计算形函数对结构坐标的导数3*8end
InterpolationMatrix=zeros(8,8);%求解节点应力应变的插值矩阵%循环高斯点for X=1:2for Y=1:2for Z=1:2E1=GaussCoordinate(X); E2=GaussCoordinate(Y); E3=GaussCoordinate(Z);GaussPointNumber = GaussPointNumber + 1;% 计算局部坐标下的形函数及形函数导数[N,NDerivative, ~] = ShapeFunction([E1 E2 E3], ElementNodeCoordinate);ElementNodeDisplacement=U(ElementNodeDOF);ElementNodeDisplacement=reshape(ElementNodeDisplacement,Dof,8);% 计算高斯点应变 GausspointStrain3_3 3*3的应变矩阵GausspointStrain3_3=ElementNodeDisplacement*NDerivative’;%3*8 8*3%把高斯点应变矩阵改写成6*1GausspointStrain=[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)]';% 计算高斯点应力GausspointStress = D*GausspointStrain-D*epsilon0;%总应变-温度应变GaussStrain(:,GaussPointNumber)=GausspointStrain;GaussStress(:,GaussPointNumber)=GausspointStress;InterpolationMatrix(K,:)=N;K=K+1;endendend%求得节点应力应变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';
上述便是热应力问题的整个求解过程，以一个两端固定约束的四棱柱为例，其在不同单位热应变的作用下的求解结果如下图所示
我的Matlab有限元编程精品课
推荐大家关注我的原创视频课程里面《Matlab有元编程从入门到精通》目前加餐到第35期。我还会持续更新，强烈推荐学习者订阅。
上新优惠价（限10名）
限时特价：449元（原价：499 元 ）
可开电子发票，赠送答疑专栏
提供vip群交流，课程可反复回看
本课程为matlab有限元编程专题课，课程主要以案例的形式进行讲解，中间会穿插案例中所涉及到的有限元基本理论，案例不局限于力学问题的有限元求解，还会涉及传热学、电学等问题的有限元求解。
