热弹性力学有限元基础理论以及代码
热弹性力学基础理论
力可以使物体发生变形，温度的变化也可以，即常说的“热胀冷缩”现象，相应的变形称为热变形。当温度变化引起的变形受到外界约束或自身约束时（梯度材料、异质结构等）将会产生应力，称为热应力。热弹性力学主要研究弹性范围内的热变形和热应力现象，并且假设变形对温度场的影响可以忽略不计。
热弹性的研究最早可以追溯到 19 世纪上半叶，J. M. C. Duhamel 于 1835 年 2 月 23 日在法国科学院朗读了第一篇关于热弹性的论文，并于 1837 年发表在《巴黎理工学院学报》上。热弹性力学的发展史可参考热弹性相关教材（见文末）。本文主要讨论线性热应力理论（Duhamel–Neumann）及其有限元方法。
物体的热膨胀行为可用线膨胀系数 来衡量，即当温度改变 1℃ 时，物体长度的变化和原温度下长度的比值，单位为 ，计算公式如下：
线性热应力理论认为温度变化引起的应变可直接与外力引起的应变线性叠加，由此可得到热应变公式如下：
式中 、 、 为热（正）应变， 、 、 为剪应变， 、 、 为热应力， 为线膨胀系数， 为温差， 为弹性模量， 为泊松比。
由上式可得到热应力表达式，详细推导过程如下：
(1)+(2)+(3)可得：
由(1)变形：
将(4)代入上式可得:
其中 为剪切模量， 为拉梅常数， 为热应力系数。
同理可得到 、 的表达式，如下所示：
剪应变公式不变。
将热应力公式代入平衡微分方程即可得到热弹性力学的平衡微分方程，原平衡微分方程如下：
将 代入(5)：
其中为 Laplace 算子。同理可得：
热弹性力学有限元格式
由于热弹性体的应变是外力引起的弹性应变和温度变化引起的热应变之和，即：
其中 为弹性矩阵， 为应变矩阵。 为初始（热）应变， 计算公式为：
由最小势能原理，参考《可重复使用运载器热防护系统热/力耦合数值计算研究》，可得单元刚度方程：
其中 和 与弹性分析里的一致，而 中多了热载荷项，计算公式如下：
其中 为单元上的力载荷， 为单元上的热载荷。由上式可发现，实际计算时可先进行热分析得到温度场，再转化为热负荷叠加到载荷列阵上，相当于只考虑温度对结构的耦合作用，该方法仅适用于机械载荷和热载荷变化不剧烈的时候。
得到单元刚度方程后组装得到整体刚度方程，求解线性方程组即可获取全局位移场。应变和应力的计算与弹性分析一致。
热力耦合问题求解的两种思路：
顺序（间接）耦合：按一定的顺序求解物理场，并将前一个物理场的求解结果作为载荷输入下一个物理场的求解过程。
完全（直接）耦合：根据耦合机理，将多个物理场的控制方程同时纳入求解过程，使用耦合单元，考虑所有物理变量的自由度，同时求解多个物理场。
代码
测试模型如下：
环境配置如下：
conda create -n couple3d python=3.8
conda activate couple3d
pip install meshio, pygmsh, pyvista, scipy
mesh.py 代码如下，主要包含与网格相关的操作，包括计算（热）弹性应变矩阵、（热）弹性矩阵、刚度矩阵、传热矩阵、弹性边界条件处理、热边界条件处理、温度载荷转换以及后处理。
import numpy as np
class tetra(object):
def __init__(self, vertices, elements, E, mu, k, alpha):
self.vertices = vertices # N × 3
self.elements = np.array(elements).astype(int) # M × 4
self.num_v = self.vertices.shape[0]
self.num_e = self.elements.shape[0]
self.E = E # 1 弹性模量
self.mu = mu # 1 泊松比
self.k = k # 1 热导率
self.alpha = alpha # 1 线膨胀系数
self.ele_V = [] # M × 1
self.B_elas = [] # M × 6 × 12
self.B_them = [] # M × 3 × 4
self.D_elas = [] # 6 × 6
self.D_them = [] # 3 × 3
self.K_elas = np.zeros((3 * self.num_v, 3 * self.num_v)) # 3N × 3N
self.K_them = np.zeros((self.num_v, self.num_v)) # N × N
self.force = np.zeros(3 * self.num_v) # 3N 力载荷
self.P = np.zeros(self.num_v) # N 热载荷
self.calc_volume()
# 计算四面体单元体积
def calc_volume(self):
for ele in self.elements:
volume = np.abs(np.linalg.det(np.hstack((self.vertices[ele], np.ones((4, 1)))))) / 6
self.ele_V.append(volume)
# 应变矩阵
def calc_B_elas(self):
if len(self.ele_V) == 0:
self.calc_volume()
for id in range(self.num_e):
A = np.hstack((np.ones((4, 1)),self.vertices[self.elements[id]]))
beta = [(-1)**(i + 1) * np.linalg.det(np.delete(np.delete(A, i, 0), 1, 1)) for i in range(4)]
gama = [(-1)**(i + 2) * np.linalg.det(np.delete(np.delete(A, i, 0), 2, 1)) for i in range(4)]
delta = [(-1)**(i + 1) * np.linalg.det(np.delete(np.delete(A, i, 0), 3, 1)) for i in range(4)]
self.B_elas.append(1. / (6. * self.ele_V[id]) * np.array(
[[beta[0], 0., 0., beta[1], 0., 0., beta[2], 0., 0., beta[3], 0., 0. ],
[0., gama[0], 0., 0., gama[1], 0., 0., gama[2], 0., 0., gama[3], 0. ],
[0., 0., delta[0], 0., 0., delta[1], 0., 0., delta[2], 0., 0., delta[3]],
[gama[0], beta[0], 0., gama[1], beta[1], 0., gama[2], beta[2], 0., gama[3], beta[3], 0. ],
[0., delta[0], gama[0], 0., delta[1], gama[1], 0., delta[2], gama[2], 0., delta[3], gama[3] ],
[delta[0], 0., beta[0], delta[1], 0., beta[1], delta[2], 0., beta[2], delta[3], 0., beta[3] ]]))
self.B_elas = np.array(self.B_elas)
# 热应变矩阵
def calc_B_them(self):
if len(self.ele_V) == 0:
self.calc_volume()
for id in range(self.num_e):
A = np.hstack((np.ones((4, 1)),self.vertices[self.elements[id]]))
beta = [(-1)**(i + 1) * np.linalg.det(np.delete(np.delete(A, i, 0), 1, 1)) for i in range(4)]
gama = [(-1)**(i + 2) * np.linalg.det(np.delete(np.delete(A, i, 0), 2, 1)) for i in range(4)]
delta = [(-1)**(i + 1) * np.linalg.det(np.delete(np.delete(A, i, 0), 3, 1)) for i in range(4)]
self.B_them.append(1. / (6. * self.ele_V[id]) * np.array([beta,gama,delta]))
self.B_them = np.array(self.B_them)
# 弹性矩阵
def calc_D_elas(self):
a = self.E / (1 + self.mu) / (1 - 2 * self.mu)
b = 1. - self.mu
c = (1. - 2 * self.mu) / 2.
self.D_elas = a * np.array([[b, self.mu, self.mu, 0., 0., 0.],
[self.mu, b, self.mu, 0., 0., 0.],
[self.mu, self.mu, b, 0., 0., 0.],
[0., 0., 0., c, 0., 0.],
[0., 0., 0., 0., c, 0.],
[0., 0., 0., 0., 0., c]])
# 热弹性矩阵
def calc_D_them(self):
self.D_them = self.k * np.eye(3)
# 计算单元刚度矩阵并组装整体刚度矩阵
def calc_K_elas(self):
self.calc_B_elas()
self.calc_D_elas()
for id in range(self.num_e):
Ke = self.ele_V[id] * np.dot(np.dot(self.B_elas[id].T, self.D_elas), self.B_elas[id]) # 单元刚度矩阵
# 整体刚度矩阵中的系数
v_index = np.array([[i*3+0, i*3+1, i*3+2] for i in self.elements[id]]).ravel()
# 组装整体刚度矩阵
self.K_elas[np.ix_(v_index,v_index)] += Ke
# 计算单元传热矩阵并组装整体传热矩阵
def calc_K_them(self):
self.calc_B_them()
self.calc_D_them()
for id in range(self.num_e):
Ke = self.ele_V[id] * np.dot(np.dot(self.B_them[id].T, self.D_them), self.B_them[id]) # 单元传热矩阵
self.K_them[np.ix_(self.elements[id],self.elements[id])] += Ke # 组装整体传热矩阵
# 弹性边界条件
def BC_elas(self, bc, bcType):
# 指定位移 化零置一法
# bc1 = np.array([[0,1,-10.],[1,1,-10.],[2,1,-10.]]) # [nodeID,dir,disp]
# dir x:0 y:1 z:2
if bcType == 'disp':
# 刚度矩阵修改
self.K_elas[3*bc[:,0].astype(int)+bc[:,1].astype(int),:] = 0.
self.K_elas[:,3*bc[:,0].astype(int)+bc[:,1].astype(int)] = 0.
self.K_elas[3*bc[:,0].astype(int)+bc[:,1].astype(int), 3*bc[:,0].astype(int)+bc[:,1].astype(int)] = 1.
# 载荷列阵修改
self.force[3*bc[:,0].astype(int)+bc[:,1].astype(int)] = bc[:,2]
# 指定集中力
# bc2 = np.array([[0,1,-10.],[1,1,-10.],[2,1,-10.]]) # [nodeID,dir,force]
# dir x:0 y:1 z:2
if bcType == 'force':
self.force[3*bc[:,0].astype(int)+bc[:,1].astype(int)] = bc[:,2]
# 热边界条件
def BC_them(self, bc, bcType):
# 第一类边界条件 指定温度 化零置一法
# bc1 = np.array([[0.0,100.0],[1.0,100.0],[2.0,100.0]]) # [nodeID,temp]
if bcType == 'Dirichlet':
# 传热矩阵修改
self.K_them[bc[:,0].astype(int),:] = 0.
self.K_them[:,bc[:,0].astype(int)] = 0.
self.K_them[bc[:,0].astype(int), bc[:,0].astype(int)] = 1.
# 载荷列阵修改
self.P[bc[:,0].astype(int)] = bc[:,1]
# 第二类边界条件 指定热流密度
# bc2 = np.array([[0.0,1.0,2.0,50.0],[1.0,2.0,4.0,50.0]]) # [nodeID1,nodeID2,nodeID3,q]
if bcType == 'Neumann':
vertices = self.vertices[bc[:,:3].astype(int)]
# 三角形面积
areas = 0.5 * np.linalg.norm(np.cross(vertices[:,1]-vertices[:,0], vertices[:,2]-vertices[:,0]),axis=1)
# 载荷列阵修改
for i in range(bc.shape[0]):
self.P[bc[i,:3].astype(int)] -= bc[i,3]*areas[i]/3*np.array([1,1,1])
# 第三类边界条件 指定环境温度和对流换热系数
# bc3 = np.array([[0.0,1.0,2.0,10.0,100.0],[1.0,2.0,4.0,10.0,100.0]]) # [nodeID1,nodeID2,nodeID3,h,T]
if bcType == 'Robin':
vertices = self.vertices[bc[:,:3].astype(int)]
# 三角形面积
areas = 0.5 * np.linalg.norm(np.cross(vertices[:,1]-vertices[:,0], vertices[:,2]-vertices[:,0]),axis=1)
for i in range(bc.shape[0]):
# 传热矩阵修改
self.K_them[np.ix_(bc[i,:3].astype(int),bc[i,:3].astype(int))] += bc[i,3]*areas[i]/12*np.array([[2,1,1],[1,2,1],[1,1,2]])
# 载荷列阵修改
self.P[bc[i,:3].astype(int)] += bc[i,3]*bc[i,4]*areas[i]/3*np.array([1,1,1])
# 温度载荷
def TempLoad(self, temp, temp0):
# 由节点温度得到温度变化以及单元温度变化
dT = temp - temp0
dT_e = dT[self.elements].mean(axis=1)
# 逐单元计算节点热载荷并加到整体载荷列阵上
for id in range(self.num_e):
strain_e = self.alpha * dT_e[id] * np.array([[1.,1.,1.,0.,0.,0.]])
dF = self.ele_V[id] * np.dot(np.dot(self.B_elas[id].T,self.D_elas),strain_e.T)
self.force[3*np.repeat(self.elements[id],repeats=3) + [0,1,2]*4] += dF.reshape(-1)
# 后处理计算应力
def post_stress(self, disp):
disp_e = np.array([[disp[3*i[0]],disp[3*i[0]+1],disp[3*i[0]+2],
disp[3*i[1]],disp[3*i[1]+1],disp[3*i[1]+2],
disp[3*i[2]],disp[3*i[2]+1],disp[3*i[2]+2],
disp[3*i[3]],disp[3*i[3]+1],disp[3*i[3]+2]] for i in self.elements])
# 单元应力
stress_e = np.array([np.dot(np.dot(self.D_elas, self.B_elas[i]), disp_e[i]) for i in range(self.num_e)])
# 单元mises应力
stress_mises_e = np.sqrt(0.5*((stress_e[:,0]-stress_e[:,1])**2 + (stress_e[:,1]-stress_e[:,2])**2 +
(stress_e[:,2]-stress_e[:,0])**2 + 6 * (stress_e[:,3]**2 + stress_e[:,4]**2 + stress_e[:,5]**2)))
# 平均法计算节点应力
stress_mises_v = np.zeros((self.num_v,2))
for i in range(self.num_e):
stress_mises_v[self.elements[i],0] += stress_mises_e[i]
stress_mises_v[self.elements[i],1] += 1
stress_mises_v = stress_mises_v[:,0] / stress_mises_v[:,1]
return stress_e,stress_mises_e,stress_mises_v
测试代码 3Dtest.py，使用 pygmsh 创建立方体模型，假设初始温度为 0℃。
import meshio
import pygmsh
import pyvista as pv
import numpy as np
from mesh import tetra
from inpTools import inp_write
import scipy.sparse as sp
import scipy.sparse.linalg as spl
# 创建模型并划分网格
with pygmsh.geo.Geometry() as geom:
# 创建立方体
geom.add_box(0, 50, 0, 50, 0, 50, mesh_size=10)
# 创建网格
mesh = geom.generate_mesh()
inp_write(mesh.points,mesh.cells_dict["tetra"],"tetrahedron.inp")
# 转为pyvista的UnstructuredGrid并可视化
vertices = mesh.points
cells = np.hstack((np.full((len(mesh.cells_dict["tetra"]), 1), 4), mesh.cells_dict["tetra"])).astype(int)
grid = pv.UnstructuredGrid(cells, [[pv.CellType.TETRA]*cells.shape[0]], vertices)
# grid = grid.explode(0.5) # 网格爆炸图
# p = pv.Plotter()
# p.add_mesh(grid,show_edges=True)
# labels = [str(i) for i in range(grid.n_points)]
# p.add_point_labels(grid.points, labels, point_size=10, font_size=10)
# p.show_axes()
# p.show()
E = 114000.0 # 弹性模量
mu = 0.3
k = 6.8
alpha = 9.2E-06
# 创建四面体网格模型
points = mesh.points
tetras = mesh.cells_dict["tetra"]
mesh = tetra(points,tetras,E,mu,k,alpha)
# 热分析 ------------------------------------------------------
mesh.calc_K_them()
# 引入边界条件
bc1 = np.array([[3.0,23.0,77.0,10.0,100.0],
[3.0,77.0,31.0,10.0,100.0],
[23.0,70.0,77.0,10.0,100.0],
[23.0,22.0,70.0,10.0,100.0],
[22.0,57.0,70.0,10.0,100.0],
[22.0,21.0,57.0,10.0,100.0],
[21.0,65.0,57.0,10.0,100.0],
[21.0,20.0,65.0,10.0,100.0],
[20.0,74.0,65.0,10.0,100.0],
[20.0,1.0,74.0,10.0,100.0],
[1.0,11.0,74.0,10.0,100.0]]) # [nodeID1,nodeID2,nodeID3(使用时需转为int),h,T]
mesh.BC_them(bc1,bcType='Robin')
bc2 = np.array([[5.0,194.0,43.0,10.0,0.0],
[5.0,48.0,194.0,10.0,0.0],
[48.0,185.0,194.0,10.0,0.0],
[48.0,49.0,185.0,10.0,0.0],
[49.0,177.0,185.0,10.0,0.0],
[49.0,50.0,177.0,10.0,0.0],
[50.0,190.0,177.0,10.0,0.0],
[50.0,51.0,190.0,10.0,0.0],
[51.0,197.0,190.0,10.0,0.0],
[51.0,7.0,197.0,10.0,0.0],
[7.0,55.0,197.0,10.0,0.0]]) # [nodeID1,nodeID2,nodeID3(使用时需转为int),h,T]
mesh.BC_them(bc2,bcType='Robin')
# 求解
temp = np.linalg.solve(mesh.K_them, mesh.P)
# # 热弹性分析 ------------------------------------------------------
# # 计算刚度矩阵
mesh.calc_K_elas()
# 热载荷施加
# 假设初始温度为 0
temp0 = np.zeros(mesh.num_v)
mesh.TempLoad(temp,temp0)
# 处理弹性边界条件
# 底面固定(z=0的所有点)
ids_bottom = np.where(points[:, 2] == 0)[0]
ids_up = np.where(points[:, 2] == 50)[0]
bc1 = np.array([[[i,0,0.],[i,1,0.],[i,2,0.]] for i in ids_bottom]).reshape(-1,3)
mesh.BC_elas(bc1,bcType='disp')
bc2 = np.array([[[i,0,0.],[i,1,0.],[i,2,0.]] for i in ids_up]).reshape(-1,3)
mesh.BC_elas(bc2,bcType='disp')
# 转为稀疏矩阵
mesh.K_elas = sp.csr_matrix(mesh.K_elas)
# 求解
disp, info = spl.cg(mesh.K_elas, mesh.force)
stress_e,stress_mises_e,stress_mises_v = mesh.post_stress(disp)
# # 创建新的UnstructuredGrid用于存储变形后模型
new_points = points.copy()
new_points[:, 0] += disp[::3]
new_points[:, 1] += disp[1::3]
new_points[:, 2] += disp[2::3]
new_grid = pv.UnstructuredGrid(cells, [[pv.CellType.TETRA]*cells.shape[0]], new_points)
new_grid.point_data['Temperature'] = temp
# new_grid.point_data['Disp X'] = disp[::3]
# new_grid.point_data['Disp Y'] = disp[1::3]
# new_grid.point_data['Disp Z'] = disp[2::3]
new_grid.point_data['Disp'] = np.sqrt(disp[::3]**2 + disp[1::3]**2 + disp[2::3]**2)
# new_grid.cell_data['Stress X'] = stress_e[:,0]
# new_grid.cell_data['Stress Y'] = stress_e[:,1]
# new_grid.cell_data['Stress Z'] = stress_e[:,2]
# new_grid.cell_data['Stress XY'] = stress_e[:,3]
# new_grid.cell_data['Stress YZ'] = stress_e[:,4]
# new_grid.cell_data['Stress ZX'] = stress_e[:,5]
# new_grid.cell_data['Von Mises Stress (Element)'] = stress_mises_e
new_grid.point_data['Von Mises Stress'] = stress_mises_v
# 使用pyvista进行可视化
p = pv.Plotter()
# 原模型
# p.add_mesh(grid,opacity=0.3,label="Original",color="grey")
# p.add_mesh(new_grid, scalars='Temperature', show_edges=True, cmap='coolwarm',
# scalar_bar_args={'vertical':True,'height':0.8,'fmt':"%.4f"},opacity=1,label="Deformed")
# p.add_mesh(new_grid, scalars='Disp', show_edges=True, cmap='coolwarm',
# scalar_bar_args={'vertical':True,'height':0.8,'fmt':"%.4f"},opacity=1,label="Deformed")
p.add_mesh(new_grid, scalars='Von Mises Stress', show_edges=True, cmap='coolwarm',
scalar_bar_args={'vertical':True,'height':0.8,'fmt':"%.2f"},opacity=1,label="Deformed")
p.show_axes()
p.show()
测试结果：
| 计算结果（Python） | 计算结果（ABAQUS） | 误差 |
|---|---|---|
|
| 2.6% |
|
| 1.1% |
） 参考文献
李维特, 黄保海, 毕仲波. 热应力理论分析及应用[M]. 北京: 中国电力出版社, 2004. 孔祥谦. 热应力有限单元法分析[M]. 上海: 上海交通大学出版社, 1999. 王洪纲. 热弹性力学概论[M]. 北京: 清华大学出版社, 1989. 马玉娥.可重复使用运载器热防护系统热/力耦合数值计算研究[D].西北工业大学,2006.
来源：有限元术
登录后免费查看全文
著作权归作者所有，欢迎分享，未经许可，不得转载
首次发布时间：2024-11-22
最近编辑：5天前
硕士
|
cae工程师
签名征集中
相关推荐
最新文章
热门文章
