有限元进阶编程——温度应力

我们之前的讨论只考虑了由外力引起的应力和变形，实际结构中，的应力和变形还会收到温度变化的影响。温度引起的应力变化我们习惯上称为温度应力。

从物理学角度分析，温度变化将会导致结构的变形，结构的变形同样也会产生热量，从而引起温度变化，两者相互交替，若全面考虑温度影响，势必将十分复杂。“幸运”的是，工程中大多数问题，变形对温度的影响很小，可以忽略不计。如此一来，就可以把温度场看做成一个独立的问题来处理，将温度应力问题看成由于温度变化导致的初应变问题。

本次推文将从理论公式入手，带着大家理解考虑温度应力后的有限元方程将会有哪些变化？最后给出相应的 Matlab 代码，加深理解。

声明：本文的理论基础及代码样例源自《结构分析的有限元方法与MATLAB程序设计》，对于学习有限元编程有着很大的帮助，课后均有源码对照练习，可在后台回复：结构有限元，即可自动获取书中全部源代码。

理论基础

加入初应变影响后的本构方程：

  初应变      可表示为：

  其中，     是材料的热膨胀系数向量，     是温度的变化量。对于各向同性材料，有

将温度初应变看成另一种荷载，利用最小势能原理，将其写成等效节点力，

温度应力问题和热荷载的应力分析问题相比，除增加一项以初应变形式出现的温度荷载以外，则完全相同。通常的处理过程是：由实验或计算得到的结构的温度场，然后把温度作为热荷载进行结构的应力分析。

对于弹性模量为      ，截面积为      ，梁长为      的梁单元，上式可以简化为

一般假设温度场沿梁长线性分布，上式可演化为

其中，    、    分别是梁单元两个节点处的温度变化量。

Matlab 程序

function etf = EquivalentThermalForce( ie )
%  计算单元的温度荷载的等效节点力
%  输入参数
%      ie  ----- 节点号
%  返回值
%      etf ----- 整体坐标系下的等效节点力 
    global gElement gNode gMaterial
    
    dT1 = gNode( gElement( ie, 1 ), 3 ) ;
    dT2 = gNode( gElement( ie, 2 ), 3 ) ;
    E = gMaterial( gElement( ie, 3 ), 1 ) ;
    A = gMaterial( gElement( ie, 3 ), 3 ) ;
    alpha = gMaterial( gElement( ie, 3 ), 5 ) ;
    Nx = E*A*alpha ;
    
    etf = [ -Nx; 0; 0; Nx; 0; 0 ] ;
    T = TransformMatrix( ie ) ;
    etf = T * etf ;
    etf = (dT1+dT2)/2 * etf ;
return

这本书的代码风格采用了全局变量形式，节点温度变化量储存在节点信息gNode第三列，在处理节点荷载时：

for ie=1:1:element_number
    egf = EquivalentGravityForce( ie ) ;
    etf = EquivalentThermalForce( ie ) ;
    i = gElement( ie, 1 ) ;
    j = gElement( ie, 2 ) ;
    f( (i-1)*3+1 : (i-1)*3+3 ) = f( (i-1)*3+1 : (i-1)*3+3 ) + egf( 1:3 ) + etf( 1:3 );
    f( (j-1)*3+1 : (j-1)*3+3 ) = f( (j-1)*3+1 : (j-1)*3+3 ) + egf( 4:6 ) + etf( 4:6 );
end

egf暂且不用管，是考虑自重影响，从以上代码片可以看出有限元方程是怎样处理温度应力（将温度初应变看成另一种形式的等效节点荷载）。

以上就是在有限元分析时，温度应力的处理方式，希望对初学有限元编程的你，有所帮助。