基于ANSYS的自适应网格划分(原创案例，转载请注明出处，谢谢！技术邻ID有限元中解人生)

基于ANSYS的自适应网格划分

何为网格自适应划分？

ANSYS程序提供了近似的技术自动估计特定分析类型中因为网格划分带来的误差。（误差估计在ANSYS Basic Analysis Procedures Guide第五章中讨论。）通过这种误差估计，程序可以确定网格是否足够细。如果不够的话，程序将自动细化网格以减少误差。这一自动估计网格划分误差并细化网格的过程就叫做自适应网格划分，然后通过一系列的求解过程使得误差低于用户指定的数值（或直到用户指定的最大求解次数）。

 

自适应网格划分的先决条件

ANSYS软件中包含一个预先写好的宏，ADAPT.MAC，完成自适应网格划分的功能。用户的模型在使用这个宏之前必须满足一些特定的条件。（在一些情况下，不满足要求的模型也可以用修正的过程完成自适应网格划分，下面还要讨论。）这些要求包括：

(1)    标准的ADAPT过程只适用于单次求解的线性静力结构分析和线性稳态热分析。

(2)    模型最好应该使用一种材料类型，因为误差计算是根据平均结点应力进行的，在不同材料过渡位置往往不能进行计算。而且单元的能量误差是受材料弹性模量影响的。因此，在两个相邻单元应力连续的情况下，其能量误差也可能由于材料特性不同而不一样。在模型中同样应该避免壳厚突变，这也可能造成在应力平均是发生问题。

(3)    模型必须使用支持误差计算的单元类型。（见表1）

(4)    模型必须是可以划分网格的：即模型中不能有引起网格划分出错的部分。

 

表1  自适应网格划分可用单元

 

2-D Structural Solids

PLANE2  2-D 6-Node Triangular Solid

PLANE25 Axisymmetric Harmonic Solid

PLANE42 2-D 4-Node Isoparametric Solid

PLANE82 2-D 8-Node Solid

PLANE83 Axisymmetric Harmonic 8-Node Solid

 

3-D Structural Solids

SOLID45 3-D 8-Node Isoparametric Solid

SOLID64 3-D Anisotropic Solid

SOLID73 3-D 8-Node Solid with Rotational DOF

SOLID92 3-D 10-Node Tetrahedral Solid

SOLID95 3-D 20-Node Isoparametric Solid

 

3-D Structural Shells

SHELL43 Plastic quadrilateral Shell

SHELL63 Elastic Quadrilateral Shell

SHELL93 8-Node Isoparametric Shell

 

2-D Thermal Solids

PLANE35 2-D 6-Node Triangular Solid

PLANE75 Axisymmetric Harmonic Solid

PLANE55 2-D 4-Node Isoparametric Solid

PLANE77 2-D 8-Node Solid

PLANE78 Axisymmetric Harmonic 8-Node Solid

 

3-D Thermal Solids

SOLID70 3-D 8-Node Isoparametric Solid

SOLID87 3-D 10-Node Tetrahedral Solid

SOLID90 3-D 20-Node Isoparametric Solid

 

3-D Thermal Shells

SHELL57 Plastic Quadrilateral Shell

 

如何使用自适应网格划分：基本过程

进行自适应网格划分的基本过程包括如下步骤：

1.       象其他线性静力分析或稳态热分析一样，先进入前处理器（/PREP7或Main Menu>Preprocessor）。然后指定单元类型，实参和材料特性，要满足上面提到的条件。

2. 用实体建模过程建立模型，用可以划分网格的面或体建模。用户不需指定单元大小也不用划分网格，ADAPT宏会自动划分网格。（如果要同时划分面和体网格，生成ADAPTMSH.MAC用户子程序－见后。）

3. 在PREP7中或在SOLUTION(/SOLU或Main Menu>Solution)中指定分析类型，分析选项，载荷和载荷步选项。在一个载荷步中仅施加实体模型荷载和惯性荷载（加速度，角加速度和角速度）。（通过ADAPTBC.MAC用户子程序       可以施加有限单元载荷，固连和约束方程。通过ADAPTSOL.MAC用户子程序       可以加入多个载荷步。这些子程序在后面还要讨论。）

4. 如果在PREP7中，退出前处理器[FINISH]。（可以在SOLUTION或在初始状态下激活ADAPT宏）。

5. 用下列方法激活自适应求解。

   Command: ADAPT

   GUI: Main Menu>Solution>Adaptive Mesh

注意，可以在热或结构分析中使用ADAPT宏，但不能在一次自适应分析中同时进行这两种不同类型的计算。在自适应网格划分的迭代过程中，单元的大小将作调整（在FACMN和FACMX指定的范围内）以减小或增加单元能量误差，直到误差满足指定的数值（或指定的最大求解次数）为止。

 6. 当自适应网格计算收敛时，程序自动将单元形状检查打开[SHPP,ON]。然后返回SOLUTION或初始状态，这取决于激活ADAPT的状态。接下来可以进       入POST1用标准操作进行后处理。

修改基本过程

 

选择自适应性

如果用户清楚某个部分网格划分的误差相对影响较小时（如应力水平较低且变化较小），可以将这些区域从自适应网格划分中排除以加快分析速度。同样，用户也许想将接近应力奇异点的部分（如集中载荷）排除掉。选择逻辑操作可以解决这类问题。

 

图 选择自适应能改进有应力集中的模型

 

如果用户选择了一个关键点集，ADAPT宏仍将包含进所有的关键点（在选择的和未选择的关键点都作网格改动），除非将ADAPT命令（Main Menu>Solution>Adaptive Mesh）中KYKPS设为1。

如果用户选择了一个面或体集，ADAPT宏将只在选择的区域调整网格大小。此时必须在激活ADAPT宏之前在PREP7中对整个模型进行网格划分。

 

用用户子程序定制ADAPT宏

标准的ADAPT宏并不能满足特定的分析需要。例如，用户可能想同时对面和体进行网格划分，这在标准宏当中是不可以的。对于这种或其他一些类似情况，可以对ADAPT宏进行修改使之适用于特定的分析。ANSYS程序用宏这种方式完成自适应网格划分，本身就使得用户可以对其进行相应的修改以适应不同的要求。方便的是，用户不用总是通过修改ADAPT代码的方式来定制宏。宏的三个部分可以用用户子程序的方法来修改，这个方法将ADAPT宏和用户文件分开，用户可以生成子程序由ADAPT宏来调用。这三个部分是：网格划分命令序列，边界条件命令序列和求解命令序列。相应的用户子程序名为ADAPTMSH.MAC,ADAPTBC.MAC和ADAPTSOL.MAC。下面看一下这三个子程序的功能：

 

生成用户网格划分子程序（ADAPTMSH.MAC）

缺省情况下，如果模型中包含一个或多个体时，ADAPT宏将只对体划分网格而不对面进行划分。如果当前选择集中没有体，宏才对面进行划分。如果要同时对面和体进行划分的话，就要生成用户子程序ADAPTMSH.MAC来提供相应的操作。在重新划分网格之前要清楚所有这些实体中划分过网格的实体。子程序大致如下：

C*** Subroutine ADATMSH.MAC-Your name-Job Name-Data Created

TYPE,1             !指定划分网格的单元类型属性

ACLEAR,3,5,2       !清除本程序中要重新划分网格的面和体的网格

VCLEAR,ALL

AMESH,3,5,2        !对面3和5划分网格（ADAPT不对其他面划分网格）

TYPE,2             !改变单元类型，划分体

VMESH,ALL          !对所有体划分网格

请查阅TYPE,ACLEAR,VCLEAR,AMESH和VMESH命令得到更详细的解释。

     我们建议将C***行包含在文件中以区别不同的宏。这一行将在任务的输出中出现，就可以确认ADAPT宏正确地调用了用户子程序。

 

生成用户边界条件子程序（ADAPTBC.MAC）

ADAPT宏在每次求解循环时都要清除并重新划分网格，因此模型的结点和单元也要不断的改变。这就使得所有的有限单元载荷，自由度固连，约束方程等施加在结点和单元的边界条件都不能使用。如果要包含这些有限单元约束时，就要使用用户子程序，ADAPTBC.MAC。在这个子程序中，可以选择结点然后定义有限单元载荷，自由度固连和约束方程。下面是一个ADAPTBC.MAC的例子：

C*** Subroutine ADAPTBC.MAC-Your name-Hob Name-Data CReated

NSEL,S,LOC,X,0      !选择X坐标为0的结点

D,ALL,UX,0          !令选择的所有结点UX为0

NSEL,S,LOC,Y,0      !选择Y坐标为0的结点

D,ALL,UX,0          !令选择的所有结点UY为0

NSEL,ALL            !选择所有结点

 

生成用户求解子程序（ADAPTSOL.MAC）

ADAPT宏中的求解命令序列很简单：

/SOLU

SOLVE

FINISH

这个缺省的命令序列只能求解单个载荷步。用户可以将其他的命令序列加入用户子程序ADAPTSOL.MAC中。

 

关于用户子程序的一些其他说明

用户可以象生成其他文件一样生成这些子程序。就是说可以用APDL语言的*CREATE命令（Utility Menu>Macro>Create Macro）和APDL的*END命令，也可以用外部的文件编辑器。当ADAPT宏调用这些子程序时，软件先搜索ANSYS根目录，再搜索用户根目录，最后是当前目录。因此，要确保其他目录中不包括与所用文件同名的文件。解释行（C***）会在输出文件中出现，可以通过它来检查是否使用了正确的文件。另外，通过在运行ADAPT宏之前用/PSEARCH,OFF(Utility Menu>Macro>Macro Search Path)指定软件只搜索ANSYS根目录和用户当前目录可以从某种程度上减少文件混用的可能。不管这些子程序在什么位置，它们都能被找到，除非将ADAPT命令的KYMAC选项设为1。

 

定制ADAPT宏(UADAPT.MAC)

有些情况下用户需要修改ADAPT宏但不能通过单独的用户子程序的方式，那么就需要直接修改ADAPT宏的主体。但是，因为某些原因，我们不推荐直接对ADAPT宏进行修改。（例如，别的用户和你同时使用一个软件，在调用ADAPT宏时会发现宏被修改了！）因此，在ANSYS安装中支持一个宏的拷贝文件UADAPT.MAC，便于用户修改。

如果对UADAPT.MAC文件进行了修改，我们建议对修改后的文件取一个新的文件名。然后在调用时输入这个文件名。要知道的是，如果新文件名是一个“unknown command”，ANSYS将搜索上级目录，然后是登录的目录，最后是工作目录，直到找到这个宏为止。如果修改的宏只能为一个用户使用，那么存储的位置应在用户登录目录的层次之下（不能等于或高于这个目录层次）。这样，存储的低层次的文件可以通过*USE命令（Utility Menu>macro>Execute Data Block）来调用。

 

自适应网格划分的一些说明

下面的建议可能有助于自适应网格划分的使用：

不需指定初始网格大小，但指定大小可能有利于自适应收敛。如果用户指定了关键点网格大小，ADAPT宏在第一次循环时使用这个值，然后在随后的循环中进行调整。用下列命令指定单元大小：

Command: KESIZE

GUI: Main Menu>Preprocessor>-Meshing-Size Cntrls>-Keypoints-All KPs

    Main Menu>Preprocessor>-Meshing-Size Cntrls>-Keypoints-Picked KPs

如果定义了线分段数或大小比例，ADAPT宏将在每次循环中都使用这个数值而不作改变。如果没有定义任何形式的网格份数，在初始网格划分时将使用缺省的网格大小[参见SMRTSIZE和DESIZE命令]。用下列方法指定线分段数或大小比例：

Command: LESIZE

GUI: Main Menu>Preprocessor>-Meshing-Size Cntrls>-Lines-All Lines

    Main Menu>Preprocessor>-Meshing-Size Cntrls>-Lines-Picked Lines

—         映射网格划分适用于2-D实体和3-D壳单元。但面的映射划分效果不明显。

—         映射网格划分适用于3-D实体。对体进行映射划分比自由划分效果要好的多。

—         总体上说，在自适应网格中有中间结点的单元比线性单元要好。

—         不要用集中载荷或尖角等引起奇异性的结构，因为此时ADAPT在这些奇异点处能量值将不收敛。如果模型中有集中载荷时，将其用施加在一个小面上的压力等效。（或通过选择将奇异部分排除在自适应网格划分之外。）

—         在许多情况下，用一系列相对小的区域替代少数几个大的区域将得到更好的网格划分。

—         如果最大响应位置已知或事先可以推测，就在附近放置一个关键点。

—         如果是在交互方式下运行ADAPT，而ANSYS在没有提示出错信息时突然退出，可以在Jobmame.ADPT文件中查看自适应网格划分部分以确定出错原因。    同样，在批处理方式下运行ADAPT时，可以看Jobname.ADPT确定出错原因。

—         如果模型中有些区域有过度的扭曲时，在网格划分中就会出错。在这种情况下，用KESIZE命令(Main Menu>Preprocessor>-Meshing-Size Cntrls>-Keypoints-Picked KPs)中SIZE域指定扭曲区域附近关键点的最大单元长度。同时，ADAPT命令中的FACMX将设为1，阻止过度扭曲部分单元大小增加。

—         应当存储结果文件（Jobname.RST或Jobname.RTH）。在ADAPT运行过程中程序如果发生中断，结果文件中将保存ADAPT过程已完成求解的内容。

—         在自适应网格运行之前应输入SAVE命令（Utility Menu>File>Save as Jobname.db）。在程序出错中断时，可以用Jobname.db重新启动计算。

《ANSYS Verification Manual》中的实例：

VM193               二维热对流的自适应分析

VM205               均布载荷下椭圆薄膜的自适应分析

下面以带圆孔矩形平板的构应力集中分析来说明自适应网格方法。

如图1所示，一个带圆孔的矩形薄板左右两边受均布拉力。

几何尺寸及材料属性如下：w=h=10mm,R=0.5mm，E=2e5MPa，μ=0.3，q=1MPa。

由于模型和载荷具有对称性，因此只需要考虑1/4模型。

下面分别是利用菜单和命令流方式进行有限元分析的方法。

图1  带圆孔矩形平板的拉伸示意图

方法一、通过加密网格实现应力集中的精确模拟

1．GUI菜单建模分析过程

第一步，清除内存准备分析。

（1）清除内存。

GUI：Utility Menu>File>Clear& Start New，单击OK按钮。

（2）更换工作文件名。

GUI：Utility Menu>File>Change Jobname，输入the plate under tension。

（3）定义标题。

GUI：Utility Menu>File>Change Title，输入the plate under tension analysis。

第二步，创建有限元模型。

（1）进入前处理器。

GUI：Main Menu>Preprocessor

（2）定义单元类型。

GUI：Main Menu>Preprocessor>Element Type>Add/Edit/Delete，选择plane42单元。

（3）定义材料。

GUI：Main Menu>Preprocessor>Material Props>Material Models，弹出对话框，连续双击Structural>linear>Elastic>Isotropic，弹出对话框，在值域EX处输入2e5，在值域PRXY处输入0.3。

（4）创建矩形。

GUI：Main Menu>Preprocessor>Modeling>Create>Areas>Rectangle，选择By Dimensions，弹出对话框，在x1、x2值域中分别输入0和5，y1、y2值域中分别输入0和5。

（5）创建实心圆。

GUI：Main Menu>Preprocessor>Modeling>Create>Areas>Circle，选择Solid Circle，弹出对话框，在WP x和WP y域中分别输入0和0，在Radius域中输入0.5。

（6）布尔运算。

GUI：Main Menu>Preprocessor>Modeling>Operate>Booleans>Subtract，选择Areas，弹出对话框，在Subtract Areas对话框中输入1，单击OK按钮，接着弹出对话框，在对话框中输入2。

（7）设置单元属性。

GUI：Main Menu>Preprocessor>Meshing>>Mesh Attributes>All Areas，弹出对话框，MAT对应项设置为1，TYPE对应项设置为1 PLANE42。

第三步，静力学分析。

（1）进入求解器。

GUI：Main Menu>Solution

（2）选择分析类型。

GUI：Main Menu>Solution-Analysis Type- New analysis，选中Static。

（3）施加对称边界条件。

GUI：Main Menu>Solution>Define Loads>Apply> Structural>Displacement>Symmetry B.C.>

On Lines，弹出对话框，选择9和10两条线。

（4）施加均布拉力。

GUI：Main Menu>Solution>Define Loads>Apply>Structural>Pressure>On Lines，弹出对话框，选择线2，然后单击OK按钮，弹出Apply PRES on lines对话框，在value栏里输入-1。

（5）自适应执行求解。

GUI：Main Menu>Solution>Solve>Adaptive Mesh，如图2所示，发现自适应网格求解选项并未出现，这时候点击Unabridged Menu如图3所示，自适应网格求解选项就出现了。点击Adaptive Mesh，弹出Adaptive Meshing and Solution对话框如图4所示，在NSOLN即迭代求解次数输入1，STARGT即结构误差输入2,其中2代表2%，FACMN即最小单元尺寸输入0.05，FACMX即最大单元尺寸输入1，点击OK按钮按钮自适应网格划分求解。

            图2  未点击Unabridged Menu之前              图3  Adaptive Mesh选项出现

图4  Adaptive Mesh自适应网格划分求解设置

（6）退出求解器。

GUI：Main Menu>Finish

第四步，执行后处理。

（1）进入后处理器。

GUI：Main Menu>Solution>General Postproc

（2）切换结果坐标系到柱坐标下。

GUI：Main Menu>Solution>General Postproc>Options for Outp，弹出对话框，选择全局柱坐标Global cylindric。

（3）柱坐标下观察结果。

GUI：Main Menu>Solution>General Postproc>Plot Results>Contour Plot-Nodal Solu，弹出对话框，选择Y-Component of stress，应力云图如图5所示。

图5  柱坐标下的σ_y分布

（4）显示单元应变能。

GUI： Main Menu>General Postproc>Plot Results>Contour Plot>Element Solu，弹出 Contour Element Solution Data对话框，选择Energy-Strain energy，应变能云图如图6所示。

图6  单元应变能云图

（5）显示单元能量偏差。

GUI： Main Menu>General Postproc>Plot Results>Contour Plot>Element Solu，弹出 Contour Element Solution Data对话框，选择Error Estimation-Structural Error Energy，单元能量偏差云图如图7所示。

                                          图7  单元能量偏差云图

（6）显示单元应力偏差。

GUI： Main Menu>General Postproc>Plot Results>Contour Plot>Element Solu，弹出 Contour Element Solution Data对话框，选择Error Estimation-Absolute Maximum Stress Variation，单元应力偏差云图如图8所示。

                                            图8  单元应力偏差云图

结果讨论：

 

                                           图9  迭代次数为1时的单元  

                                    图10  迭代次数为2时的单元

                                   图11  迭代次数为3时的单元

                                   图12  迭代次数为4时的单元

                                    图13  迭代次数为5时的单元

图14  迭代次数为6时的单元

表1 自适应网格不同迭代次数结果比较

迭代次数	节点数	单元数	symax	symin
1	49	34	2.827	-0.865265
2	394	347	3.094	-1.06
3	361	314	3.074	-1.062
4	365	318	3.087	-1.064
5	360	313	3.102	-1.063
6	354	307	3.101	-1.065

表2各种方法结果比较

方法	一般方法(esize=0.1)	一般方法(esize=0.05)	局部加密法	p单元法	自适应网格法
节点数	2418	8894	944	542	354
单元数	2318	8695	873	163	307
symax	2.975	3.077	3.107	3.092	3.101
symin	-.966125	-1.046	-1.063	-1.069	-1.065

从表1可以看出随着迭代次数的增加结果趋于稳定，应力集中系数稳定在3.10左右，单元数及节点数随着迭代次数的增加呈现波动趋势，原因在于设定的结构误差目标1%没办法满足，这点可以从output window看出。针对这种情况可以通过减小最小单元最小尺寸来实现。

从表2可以看出p单元法节点数最少，自适应网格法单元数最小，它们均在代价最小的情况下得到了非常好的结果；一般网格加密的方法必须要整体加密才能得到较好的结果，原因在于应力梯度不大的区域加密网格效果也不明显，只有在应力梯度较大的区域加密网格才能达到较好的效果，因此局部加密网格法应运而生，在应力梯度较大的位置加密网格，既能达到分析的精度又能避免浪费计算资源，可以说是针对应力集中不错的一种方法。当然局部加密法的前提是要能知道哪些区域应力较大，因此在实际操作中可以先采用一般网格法，得到一个粗糙的结果，根据结果再来有目的性的局部加密网格，从而达到较好的分析效果。

综上，在实际分析应力集中的时候可以采用p单元法、自适应网格法和局部网格加密法，这里需要说明的是针对三维实体单元p单元法以及自适应网格法计算规模可能会比较大，实体单元采用局部加密计算规模相对比较容易控制。

附件为自适应网格对于的命令流以及完整的文档，象征性收费1元。

视频链接：http://www.jishulink.com/college/video/c10819，有兴趣的可以观看，谢谢