CAE前处理 | 选择合适的网格密度：圆角
前文回顾：
01 前言
前文我们对结构刚度计算问题需要的实体单元数量进行了一个基本的探讨，通过实例计算和对比大家会发现：高阶单元和低阶六面体单元，只要按照常规去划分网格，单零件的静刚度问题计算并不困难（薄壁问题需要特别注意）。
但好像平时分析大家总是担心自己的网格精度不够，并且强调网格无关性的重要性，这又是为什么呢？
实际上，通常我们担心的并不是网格对于刚度（变形）计算精度不够，而是担心网格对于强度（应力）计算精度不够，这是由于刚度问题是一个全局问题，而强度问题是一个局部问题，这也就意味着，需要在局部划分足够的网格才能捕捉到详细的应力变化情况（关于局部加密技术见文章CAE前处理 | 网格局部加密）。
然而具体要划分多少网格才能够捕捉到局部应力呢？说实话，这个问题非常复杂并且没有定论，因为任何一个结构局部应力的变化和结构的形状以及受载方式息息相关，我们基本只能通过网格无关性验证去得知当前的网格密度是否合适。
可是，大家会发现网格无关性验证工作量巨大并且对于装配体分析往往不太现实。因此，对于一些典型的结构特征如果能够得到一个大致的参考规律，那么很多时候就不必再进行网格无关性验证或者只用进行少数的验证。
针对这样一种需求，本文包括以后一些文章会间断地针对一些典型的特征进行分析对比，意图得到哪怕一丢丢对于应力计算有用的经验结果。由于圆角是最为常用的一类几何特征，因此本文主要就圆角的应力计算进行一定的对比说明并试图得到些许规律。当然，我并不希望大家直接翻到最后看结果，而希望能够关注于探讨的过程和方法，因为就像前文所说，应力问题千千万，知道这个问题的结果并不会对自身的能力有多大帮助。
02 模型选取
与前篇文章思路一致，要得到规律性结果的前提是找到具有典型特征的结构模型和工况。经过筛选，本文选取有圆角过渡的台肩模型作为标准对比模型：
主要原因有几点：
①模型经典，在一般的应力集中系数手册中都可以查到其理论/试验结果 （这里查询大量资料貌似有限元计算结果和光测弹性力学试验得到的结果存在大概10%的差距，因此本文暂时不以试验结果做为参考标准）
②工况经典，弯曲/拉伸/扭转都是圆轴可能受到的典型工况，因此也具有典型的应力集中现象
③拓展性强，圆轴模型作为典型的广义平面应变模型，得到的很多结果可以直接拓展到2D问题中，并且由于具有3D特征，因此很多结论同样也适用于其它三维模型。
因此，本文选择该模型作为对比的标准模型。并且由于低阶四面体，高阶四面体，低阶六面体，高阶六面体都是现在常用的实体单元类型，所以文章会对这些单元的计算结果都进行一个评估对比。
03 分析规划
分析对比一定不是盲目地，而是需要针对所分析的问题去仔细规划，而正式规划之前，我们还必须对所分析的问题具有一定程度的了解。
台肩在受到标准拉伸载荷时的应力分布：
首先，应力热点位置出现在台肩过渡圆角根部，也就是说这里是后续进行局部加密的重点区域。
其次，应力集中具有轴对称性，因此后续并不需要对周向网格的划分进行探讨。
下面我们来看下结构的应力在绕轴向和径向的分布特点：
显然，绕轴线方向和径向方向上应力均具有典型的非线性变化，也就是说和常规的弯曲/扭转线性应力，拉伸均匀应力不同，我们确确实实需要使用更密的网格才能去捕捉到这些非线性变化趋势。
同时我们也会发现，沿着轴线方向的应力变化比径向方向更加附加和剧烈，所以，为了减少前处理的参变量，主要以研究轴线方向上的网格划分为主即可。
台肩在受到标准扭转载荷时的应力分布：
其实大家会发现，圆轴受到扭转载荷的应力分布比拉伸载荷要更加有规律性，同样是轴对称规律，并且轴扭转的切应力理论上规律性会更加接近线性，同样提取结构的应力在绕轴向和径向的分布特点：
和预估的一样，受到扭矩作用时结构的应力变化趋势相对更加缓和，也就是说对于台肩模型，如果能够达到拉伸载荷局部应力捕捉要求应该也是能够满足扭转载荷的需求。
台肩在受到弯曲为主载荷时的应力分布：
相比于其它两种工况，受到弯曲载荷时结构的应力分布已经不具备轴对称性质，应力分布呈现"斑状"，而且预计扭转时候应力的非线性变化程度应该是最大的，提取结构的应力在绕轴向，径向以及周向的分布特点：
可以看到，三个方向的应力变化中非线性程度：轴向＞径向＞周向，并且轴向和深度方向上的应力变化趋势以及剧烈程度和受到拉伸载荷比较接近。
现在，我们将三种受力状态下的应力分布进行一个简单的对比会发现：对于圆角特征的局部应力，变化最为剧烈的为沿着圆角方向（轴向），其次为下覆深度方向（径向），最后为环绕圆角方向（周向）。其中绕圆角周向就算按照弯曲应力工况去看，其非线性程度并不大，只要能够合理的离散几何特征即可。
虽然实际结构受力是综合受力状况，但是我们已经将典型非线性变化趋势提取出来了，因此只需要着重研究圆角方向网格划分，并以深度方向网格划分为辅即可。
这样我们的分析规划如下：
选取圆轴受到拉伸载荷为典型对比工况（具有代表性），分别使用1~n层不同类型网格进行对比，着重得到沿圆角变化方向网格划分层数规律，并且留意下覆深度方向对结构应力的影响。
04 深度方向影响
本文并不去详细对比深度方向网格的规律性变化，主要有几点原因：
①深度方向在结构内部不便于控制，就算得到了规律也很难在工程应用中使用。
②深度方向应力的变化剧烈程度不绕圆角变化方向，只要我们找到圆角方向的网格规律，按照正常的网格过渡也是完全可以满足深度方向要求。
但是，并不是说深度方向不重要，为了强调这个问题，这里对比两组模型：
这两组有限元模型除了在局部深度方向上网格密度不同，其余网格划分规律，数量以及工况载荷都完全一致，现在提取整体的应力热点以及深度方向应力变化趋势：
大家可以看到，对于深度方向没有加密的模型，局部应力热点值为1.634MPa，而深度方向加密了的模型，局部应力热点值为1.869MPa，计算下来，如果不进行加密应力峰值偏小12.6%。
如果仔细对比深度方向的变化趋势可以看到，没有加密的模型完全没有捕捉到局部的应力变化，并且这种不精确的捕捉以及会影响到最大应力的数值，因此我们更应该注意：网格局部加密过渡一定要缓！
05 圆角方向影响
现在我们进入到文章的核心部分：圆角部分划分多少层网格合适？
为了便于对比，六面体网格和四面体网格划分需要采取不同的策略：
这里不管是六面体还是四面体网格，变量均为圆角方向网格层数，对比模型五组：1，2，4，8，12层网格，对应夹角为：180°，45°，22.5°，11.25°，7.5°，层数和度数转换示意如下：
按照这样的规律，对于不同单元类型以及划分层数可以得到以下结果：
低阶六面体：
低阶四面体：
高阶四面体：
【说明：这里由于部分原因不对高阶六面体进行分析，不过由于高阶六面体精度＞高阶四面体，因此高阶四面体得到的结果同样适用于高阶六面体】
将上述计算结果列表绘制成曲线如下：
根据计算数据可以得到，对于圆角特征，如果按照5%精度为可接受精度，那么低阶六面体至少划分12层（7.5°），低阶四面体至少划分.....（不做解释），高阶四面体至少划分4层（22.5°）。
06 估测方法
虽然通过大量数据对比能够直接得到一些指导网格划分的参数，但是会遇到两个问题：
①每一种不同特征的划分参数都需要对大量的数据进行分析和对比才能得到满意的结果，这在工程上其实并不实用。
②实际模型要复杂得多，有些时候很难通过特征去得到一些规律性可以借鉴的结果。
因此，个人这里推荐两种相对更加普适的估测方法，这些方法的根源都来源于后处理部分。
方法一：应力分布法
在前面的对比中大家不难发现，合理的应力结果意味着热点周边应力的非线性分布趋势也更加光滑，比如提取高阶四面体在网格层数为4层，8层和12层时候的圆向应力变化结果：
大家明显可以看到当层数为4层时，峰值应力曲线明显非常的粗糙，到8层的时候才略微有所改善（实际8层精度也可以接受），到12层的时候已经比较光顺。因此，这里说的应力分布法的基本思想是：当热点应力周围的应力分布趋势足够光滑，那么理论上热点应力精度也较高（前提深度方向网格足够）。
方法二：热点覆盖法
当然大家会发现应力分布法还是有些麻烦的，而且怎么样算光滑不同人的评判标准并不一致，因此有了下面这种方法：热点部分（红色 区域）至少覆盖xx层网格。
这种方法第一次听说是在教主（张晔）的ansys workbench视频中，这套视频我是非常建议分析/结构工程师好好学习一番，大家可以在机械人读书笔记公 众 号中免费学习。
这里我就这种方法再细致的和大家解释一下：
我们知道应力集中局部的应力变化是非线性的，而有限元要捕捉到这样的非线性变化需要划分足够的网格数量，那么就要求局部的应力梯度变化不能太剧烈。
大家想一想，热点（红色）区域代表什么？实际代表1/n的全局应力，比如通常云图条默认是10份，如果最大应力100MPa，那就意味着90~100MPa的结果都是红色的，从局部来看这就是限制应力的梯度变化。
那么在1/10的梯度变化内，使用2~3层线性/二次单元去拟合，一般情况下是完全足够的，或者说至少精度不会太差。
比如下面是12层低阶六面体和4层高阶四面体的结果云图：
可以看到，使用低阶六面体热点区域至少要覆盖3层，使用高阶单元热点区域完美覆盖1层单元足以。
