第七讲 耦合问题有限元分析 元计算技术部

在实际生活中，我们要解决的许多问题是很多个物理场(诸如温度场，应力场，流场等)的叠加问题，而且这些物理场之间是相互影响的，比如炼钢的时候温度的高低对于应力分布会有影响，这种多个物理场相互叠加的问题就称作多场耦合问题，ELAB1.0是基于单场的物理微分方程出发，可以比较好的实现不同物理场的耦合。本讲针对实际中比较普遍的热固耦合问题进行分析计算，以下将从该类问题的物理方程，有限元分析以及具体实例的ELAB1.0实现几个方面进行介绍。 基本方程 ELAB模型向导实现 有限元脚本文件分析

热固耦合的基本方程 热固耦合问题的基本方程包括热传导问题的基本方程和固体力学的基本方程以及两者之间的耦合关系，以下以二维问题为例。 热传导基本方程： 对于线弹性小变形问题来说，固体的变形对温度的影响比较微小，可以忽略不计，因此热固耦合中热传导过程的基本方程不变： 其中，k为热传导系数 ，q为热源 线弹性固体力学基本方程： 对于线弹性问题，其本构方程将受温度变化的影响，下面给出线弹性问题的平衡方程、几何方程以及与受温度场变化影响的本构关系 :

固体场平衡方程： 固体场几何方程： 本构关系： 其中，a为热膨胀系数 ，E表示弹性模量，ν表示泊松比

传热问题的边界条件有三类： 第一类边界条件： 第二类边界条件： 第三类边界条件: 其中，q0是边界上热流的给定值 ，nx、ny、nz分别为边界表面外法线在 x、y、z方向的的方向余弦，h表示物体与周围介质的热交换系数， T0表示环境温度 。 线弹性问题有两类边界条件： 固定位移边界条件： 边界均布力载荷条件： 其中， u0表示x方向的位移， v0表示y方向的位移，T0表示x方向的边界载荷，T1表示y方向的边界载荷 针对以上理论分析，以下用ELAB1.0公式库实现的方式求解一个相应的实际算例。

热固耦合有限元分析 工程背景 平板长1米，宽0.5米，左端温度为0℃，右端温度为100℃，下端完全固定。如下图所示，求在此条件下的板的温度分布、变形和应力。板的膨胀系数为1.0e-5/℃，弹性模量为1000MPa，泊松比为0.3，热传导系数为10W/m/℃。不计板的体力和内热源。 几何模型

热固耦合ELAB1.0软件实现 工程建模 1、点击“工程向导”进入公式库 2、选择“多物理场耦合”→“热固耦合” 3、选择“坐标系”

4、选择“单元类型” 5、选择“问题类型” 6、定义工程名和工程路径，完成工程设置

b场体单元材料参数图 b场边界单元材料参数 定义材料参数 点击工具栏“参数设置”→“材料参数”，如下图所示： 材料参数对话框中设定相应的材料参数，如下图所示： a场体单元材料参数图 a场边界单元材料参数 b场体单元材料参数图 b场边界单元材料参数 c场体单元材料参数图

前处理 点击工具栏中“前处理”按钮进入GID。 注：进入GID后要进行ELAB1.0的数据转化data→problemtype→ELAB 几何建模： 首先建立一个小的矩形面，利用gid中copy命令中的拉伸功能建立如下图所示的几何模型，详细步骤可以参考《有限元分析基础与应用》相关章节。 有限元模型 施加材料属性： 在condition窗口中为a场(位移场)、b场(温度场)和c场(应力场)分别施加材料属性和边界条件，该模型只有一种材料，材料赋值如下图所示：

a场材料设置 b场材料设置 c场材料设置 施加边界条件： 温度场边界设置 位移场边界设置

工程求解 后处理 划分网格： 点击工具栏中“求解计算”按钮，完成模型的求解计算。 网格划分(网格尺寸0.04) 温度分布云图 x方向位移分布云图

应力场dxx分布云图 应力场dyy分布云图 变形云图

有限元语言描述文件 为生成该问题有限元计算的所有程序源代码，针对之前的ELAB1.0有限元分析得到的微分方程弱形式，ELAB1.0软件提供简洁的有限元语言描述文件，包括微分方程描述文件、多物理场描述文件以及求解命令流控制文件。 针对该问题的有限元描述文件包括heatxy.fde(温度场fde文件)， delxy.fde(位移场方程描述文件)，selxy.fde(应力场方程描述文件)，couple.mdi，couple.gcn 微分方程描述文件heatxy.fde(温度场fde文件) 在heatxy.fde给出单元的待求未知量，涉及到的材料参数，单元的形函数表达式，刚度矩阵表达式和载荷表达式，以及为描述刚度矩阵和载荷向量而自定义的函数。热固耦合中热传导过程的基本方程不变，因此对应的有限元文件也不变，可参考《第六讲热传导过程的有限元分析》，详细的解析见《有限元分析基础和应用》中相关章节。

未知变量对应微分方程弱形式中的变量(本构中) 微分方程描述文件delxy.fde(位移场fde文件) 在位移场方程描述文件delxy.fde中， 给出单元的待求未知量，涉及到的材料参数，单元的形函数表达式，受温度影响的刚度矩阵表达式和载荷表达式，以及为描述刚度矩阵和载荷向量而自定义的函数。 微分方程弱形式： 未知变量： DISP u v w 未知变量对应微分方程弱形式中的变量(本构中) u v w

单元刚度矩阵对应微分方程弱形式中的左端项 耦合信息： COEF tn 耦合变量对应微分方程弱形式中的变量 T 材料参数： MATE pe pv alfa fx fy rou alpha 材 料参数行对应微分方程弱形式中的变量 E ν a fx fy 单元刚度矩阵： dist = [ev_i;ev_j]*sm_i_j*fact+[ep_i;ep_i]*shear*fact 单元刚度矩阵对应微分方程弱形式中的左端项

单元载荷向量对应微分方程弱形式中的右端的第二项和第一项 单元载荷向量： load=[u_i]*f_i*vol+[ev_i]*fte_i*vol 单元载荷向量对应微分方程弱形式中的右端的第二项和第一项 多物理场描述文件couple.mdi 2dxy #a 0 2 u v fde delxy q2 #b 0 1 u fde heatxy q2 #c 0 3 dxx dyy dxy fde selxy q2 # 坐标系(二维直角坐标系) a场0个初值2个自由度 方程描述文件+单元类型和积分方法 b场0个初值1个自由度 b场0个初值3个自由度 结束标志

求解命令流控制文件couple.gcn DEFI a ell b b ell c str a b START b START a SOLVC b SOLVC a SOLVSTR c a gidres(coor0); a场采用ell算法耦合b场数据 b场采用ell算法 c场采用str算法耦合a场和b场数据 （空一行） 初始化b场 初始化a场 求解b场 求解a场 最小二乘求解c场 输出GID格式的结果数据

THANKS