我试图求解一个位移场，我知道我的离散域上的热应力，该离散域由六角单元组成。

我的第一个问题是：除了积分线性弹性方程之外，还有更简单的方法来解决位移：

$$\nabla \cdot \sigma = 0$$ 在哪里$\sigma$是应力（包含机械应力和热应力）。

我目前正在使用有限体积法求解上述方程，积分形式为

$$\int \int (\sigma \cdot \hat{n}) dA = 0 \ \ \ \ \text{ applied divergence theorem}\\ \int \int (\sigma \cdot \hat{n}) dA \approx \sum_{f=1}^N \sigma_fA_f\hat{n}_f$$ 其中积分是使用中点规则近似的，该规则读取为单元所有面上的面应力、面面积和面法线的乘积之和。面应力是通过取相邻单元应力的平均值来计算的。

应力使用各向异性材料的 duhamel-neumann 关系（请注意，我目前正在对简单的各向同性材料测试代码） 其中是使用计算的机械应变线性应变位移关系，是热应变，使用 计算 ，其中是 CTE。

$$\sigma = C:(\epsilon_m - \epsilon_t)$$ $\epsilon_m$$\epsilon_t$ $$\epsilon_t = \alpha (T-T_0)$$ $\alpha$

我可以将上述方程组表示为线性系统，其中是要求解的位移。是强迫，基本上所有的热应力项都包含在内（热应力对没有影响）。 $Au=b$$u$$b$$b$$A$

好消息是，如果我将相同的热应力分配给整个域（即，所有单元的温度相同），那么我得到的解决方案是“正确的”。正确地说，我的意思是它与商业软件 (ANSYS) 匹配。所以有些事情是正确的。

但是，如果我在整个域中的温度分布不一样，那么我得到的解决方案不仅与 ANSYS 不匹配，而且也不符合我的直觉（温度变化较大的某些单元膨胀较少甚至收缩（ !)，如计算的位移场所示，与温度变化较小或没有变化的零件相比。

我检查了我的代码，似乎我正确地实现了一切，所以我想知道我的数字是否不正确。但是，如果将整个字段初始化为相同的温度，我将得到正确的解决方案，因此某些事情必须是正确的。例如，添加热应力项是否会以某种方式影响我应该如何进行离散化？我认为不会。这个词很简单，实际上只是一个强迫词。

如前所述，热应力以数学方式体现在力矢量中，对没有影响。所以改变温度只会影响我的力矢量的大小，而我保持不变。所以这告诉我，我$b$$A$$A$$A$矩阵的公式是正确的，或者它可能公式不正确，但不会影响我对其他问题的解决方案。我认为后者不太可能，因为我还在悬臂和各种其他弯曲/压缩问题（没有热应力）上测试了我的代码，并且我得到了与分析解决方案或商业软件非常匹配的解决方案。我还进行了一些制造解决方案研究以查看收敛性并且测试成功。

正是这个热应力项引起了我的问题，而且只有当离散域中的温度不均匀时。

有什么建议？如果需要，我很乐意提供更多信息。

编辑#2：这是我现在正在使用的特定几何图形：

长度在 x 方向。我将的均匀温度应用于所有电池。下图是y方向的位移轮廓，是上图中更细化的方向。 $1900K$

+y 和 -y 表面的轮廓的最大 y 位移约为 0.0006m 和 -0.0006m。这种材料是钢，使用 CTE = 。块的尺寸在 x、y、z 方向上约为 0.6098m x 0.0047625m x 0.00508m（抱歉有多个小数。原始网格以英寸为单位创建）。3 个方向上的网格步长是恒定的，为 0.03048m（1.2 英寸）x 0.004318m（0.17 英寸）x 0.0254m（1 英寸）。$1.2\times 10^{-5} K^{-1}$

所以在这里我得到了一个与 ANSYS 给我的解决方案相匹配的解决方案。

接下来，我继续运行不同的模拟，其中我只在顶部两层的域中指定了 1900K 的温度（请参见下面白框中的区域），而域的其余部分为 298K。

在这个只有部分域经历温度变化的模拟中，我得到的解决方案与预期或 ANSYS 给我的解决方案相去甚远。

编辑#3：为上述测试用例添加位移轮廓。

这是第二种情况下 x、y 和 z 位移的 3 个位移轮廓。我的直觉告诉我，x 和 z 位移轮廓显示了预期的趋势（不确定幅度是否正确，但至少趋势是合理的）。但我不希望 y 位移的行为方式与 y 位移轮廓中显示的方式相同。