FEtch 应用 —— 热固耦合
关键词
： 瞬态 线性 热传导 弹性力学 向后差分 热固耦合 解耦合
结构内温度场发生变化时，若受到外部约束或温度场不均匀时，会产生一定的应力，称温度应力。温度应力的出现，已引起工程上普遍关注。在土木工程中，钢筋混凝土构筑物由于受到环境温度变化的影响，表面和内部会产生变形，若遇约束，会引起温度应力， 当应力达到一定值时，结构内部产生微观裂缝，甚至发展为裂缝；大体积混凝土地基和大坝结构由于水泥浇注期内水化热的作用，冷却收缩时温度应力若超过材料抗拉强度，也会呈现裂缝；道路和桥梁由于气候多变，在路面荷载和不均匀温度场共同作用下也有可能开裂受损。由此可见，不均匀温度场和温度应力计算具有重要意义，其结果可以直接为工程设计提供依据。
本节以二维弹性体的热固耦合问题为例，介绍 FEtch 系统在求解热固耦合问题中的应用。
控制方程
瞬态热传导
瞬态热传导
二维直角坐标系下，瞬态热传导服从如下偏微分方程
\rho c \frac{\partial T}{\partial t}-\frac{\partial}{\partial x}\left(k \frac{\partial T}{\partial x}\right)-\frac{\partial}{\partial y}\left(k \frac{\partial T}{\partial y}\right)=0 \ \left(\text{in } \Omega\right) \\ 其中， \rho 表示材料密度， c 为材料的比热， T 表示温度， k 是热传导系数， x 和 y 为二维直角坐标系下的坐标变量。
热弹性力学
二维直角坐标系下，忽略惯性力和体力的作用，热弹性力学方程可写成如下形式
\frac{\partial \sigma_{x x}}{\partial x}+\frac{\partial \sigma_{x y}}{\partial y} = 0 \ \left(\text{in } \Omega\right) \\
\frac{\partial \sigma_{x y}}{\partial x}+\frac{\partial \sigma_{y y}}{\partial y} = 0 \ \left(\text{in } \Omega\right) \\ 几何方程为
\varepsilon_{xx}=\frac{\partial u}{\partial x}, \quad \varepsilon_{yy}=\frac{\partial v}{\partial y}, \quad \varepsilon_{x y}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x} \\
本构方程（平面应力）为
\left(\begin{array}{c}\sigma_{x x} \\ \sigma_{y y} \\ \sigma_{x y} \end{array}\right)=\frac{E}{(1+\nu)(1-\nu)}\left(\begin{array}{ccc} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & (1-\nu) / 2 \end{array}\right) \left(\begin{array}{l}\varepsilon_{x x} \\ \varepsilon_{y y} \\ \varepsilon_{x y}\end{array}\right) -\frac{E}{(1-\nu)}\alpha T\left(\begin{array}{l} 1\\ 1 \\ 0 \end{array}\right) \\
其中， u 、 v 为位移 \boldsymbol{u} 的分量。 \varepsilon_{x x} 、 \varepsilon_{y y} 、 \varepsilon_{x y} 为应变 \boldsymbol{\epsilon} 的分量， \sigma_{x x} 、 \sigma_{y y} 、 \sigma_{x y} 为应力 \boldsymbol{\sigma} 的分量。参数 E 为弹性模量， \nu 为泊松比， \alpha 为线膨胀系数。
算例
问题介绍
边长为 0.4\ \mathrm{m}\times0.3\ \mathrm{m} 的长方形薄板，上下边的位移竖向约束，横向自由，两侧边自由。初始温度为 0 { }^{\circ}\mathrm{C} ，底部突然施加 10 { }^{\circ}\mathrm{C} 的固定温度，上边和两侧绝热，求 100\ \mathrm{s} 内薄板的温度、变形和应力的变化情况。具体参数值为 \rho=5000\ \mathrm{kg}/\mathrm{m}^3 ， c=10\ \mathrm{J}/\mathrm{kg}/^{\circ}\mathrm{C} ， k=100\ \mathrm{W}/\mathrm{m}/{ }^{\circ}\mathrm{C} ， E＝1.0×10^{10}\ \mathrm{N/m^2} ， \nu = 0.3 ， \alpha=1.0\times10^{-5}\ /^{\circ}\mathrm{C} 。
网格剖分
根据对称性，取薄板的右半边进行建模。根据所开发程序的计算要求，网格剖分采用 4 节点等参单元。共使用了 150 个单元，176 个节点。
计算结果
对计算结果稍加整理，效果如下。
温度的演化过程
随着传热的进行，温度自下而上逐渐增大。薄板从下到上逐渐膨胀。为了更好地显示变形效果，动画中的位移放大了 800 倍。
位移的演化过程
水平位移
竖向位移
薄板从下到上逐渐膨胀。最终，水平向位移成线性分布，而竖向由于受到约束，不发生位移变化。
应力的演化过程
x 方向正应力
y 方向正应力
切应力
进一步分析可以发现数值解与解析解符合得很好，进而证明算法和程序的有效性。
了解更多最新信息，欢迎关注我们的 B 站主页 ：
阅读详细文档，请访问 FEtch 系统的 技术网站 ：
Tags - FEtch 有限元自动生成系统
有任何疑问或建议，欢迎加Q群 " FEtch有限元开发系统(519166061) " 留言讨论。
点击链接加入群聊【FEtch有限元开发系统】： https:// jq.qq.com/? _wv=1027&k=eLERtaMm
大漠胡杨工作室长期开展 FEtch 系统的免费试用活动，感兴趣的朋友可以关注微信公众号" 有限元语言与编程 "，留言" 试用 "，免费获取 许可证文件 。
