首页/文章/ 详情

算得准比能算更重要!ANSYS workbench中单元阶数和积分水平对结构模拟精度的影响分析!

6月前浏览1132
本例通过一个悬臂梁的静力结构分析案例演示ANSYS workbench中单元阶数(线性或二次)、积分水平(完全积分或减缩积分)对结构模拟精度的影响。

结构梁由于长度很大,通常通过梁单元进行模拟,但本例为了说明实体单元的效果,将梁按照实体单元建模。

本文案例内容参考了庄茁老师的《abaqus非线线有限元分析与实例》一书。分析软件虽然不同,但从结果看,ANSYS的计算结果与abaqus结果基本一致。


—01—

理论计算


悬臂梁长度150mm,宽2.5mm,高5mm。梁一段固定,另一端施加5N的集中力。如下图所示:
 

梁的材料属性如下图设置所示:

在载荷P作用下,梁自由端静力分析的挠度为:
 
其中:
I=bh³/12,L为梁长度,b为梁截面宽度,h为梁截面高度。

当P=5N时,理论上计算梁自由端位移δ、固定约束端弯矩M、端部最大正应力σ和等效应力如下:

① 梁自由端的挠度:
δ=5*150³/(3*70000*2.5*5³/12)=3.086mm
② 梁固定端的弯矩:
M=PL=5*150=750N.mm
③ 梁端部的最大正应力:
σ=M/Wz=M/(bh²/6)=750/(2.5*5*5/6)=72MPa
④ 最大等效应力:
端部剪应力很小,可按剪力除以面积近似计算剪应力:τ=5/(2.5*5)=0.4MPa,由此可计算最大等效应力为:
 
σₑ=(0.5*((72-0.4)^2+0.4^2+72^2))^0.5=71.8MPa


—02—

ANSYS中有限元计算


本例中,有限元计算时我们考虑单元的阶数和单元积分水平的影响。

有限元分析中,单元的积分分为完全积分和减缩积分。完全积分是指单元具有规则形状时,全部Gauss积分点的数目足以对单元刚度矩阵中的多项式进行精确积分。减缩积分顾名思义,则是减少了积分点数。

在ANSYS中,一阶六面体单元的完全积分的积分点数(SOLID185)通常为2×2×2=8个积分点,减缩积分则仅使用1个积分点。对于二阶六面体单元(SOLID186),完全积分的积分点数通常为3×3×3=27个积分点,但ANSYS中实际采用优化后的14个积分点方案,减缩积分则采用2×2×2=8个积分点。
 
完全积分单元

 
减缩积分单元

在ANSYS workbench中单元积分设置在几何部件下,首先需要在geometry下将element control设置为“manual”,然后再具体的几何部件的详细设置面板设置“full”或“reduced”,如下图所示:


有限元分析中,单元的阶数是指单元内部用于近似描述物理量(如位移、温度、压力等)分布的形函数的阶次。通常单元阶数采用线性(1阶,即Linear)或二次多项式(2阶,Quadratic)。当然也可采用更高阶次,但计算成本会很大。ANSYS中单元阶次的设置在网格划分处,如下图:

本例中,为体现不同网格密度对结果的影响对截面高度上的网格层数和梁长度方向的网格数量进行控制,并设置为参数化输入项。设置如下:



网格控制时,分为三种情况(后续通过参数化设置):
① 截面高度上2层网格
单元尺寸:2.5mm×25mm,数量:2×6,长宽比:25/2.5=10
② 截面高度上4层网格
单元尺寸:1.25mm×12.5mm,数量:4×12,长宽比:12.5/1.25=10
③ 截面高度上8层网格
单元尺寸:0.625mm×6mm,数量:8×25,长宽比:6/0.625≈10  

案例约束边界条件和载荷设置如下图:

本例为分析单元阶数和积分水平对结构模拟精度的影响,考虑的结果包括:变形、最大应力和支撑反力。为此对这三类结果进行参数化设置,作为参数化输出,如下各图所示:


① 变形输出设置:


② 等效应力输出设置:

③ 弯矩反力输出设置:

以上设置完成并计算后,回到workbench主界面,可以看到已经生成了参数化。如下图:

双击“Parameter Set”进入参数化设置,按照前述网格划分规划,设置DP0、DP1、DP2三种情况的网格尺寸,如下图:

设置完成后,更新所有设计点。即可得到某一种单元设置的结果。计算完成后,改变单元阶数和积分类型,得到另外的计算结果。各结果汇总如下:

①单元阶次为1阶,完全积分

② 单元阶次为2阶,完全积分

③ 单元阶次为1阶,减缩积分

④ 单元阶次为2阶,减缩积分


—03—

理论结果与有限元结果对比分析


3.1 位移结果对比
理论计算表明,本例悬臂梁的自由端的最大位移为3.086mm。下表是ANSYS关于本例的实体单元不同阶次和积分水平的计算结果:

表1 不同单元阶次和积分水平的位移结果

① 完全积分的线性单元
由上表可知,完全积分的线性单元得到的结果相当差,基本不可用。网格越粗糙,结果精度越差。即使网格很细(8×25),得到的自由端位移也只有理论的1.8108/3.086=58.7%.

可以看到,完全积分的线性单元,梁截面高度上的单元数量并未对结果精度的提升有较大帮助。这是因为自由端的挠度误差由“剪切自锁”造成。这是存在于所有完全积分的线性实体单元的问题。
 

剪切自锁的原因在于单元的边无法弯曲,导致单元在弯曲时过于刚硬,变形表现为“伪剪切变形”。其原理解释详见《有限元仿真中的两大“锁死”难题:剪切自锁与体积自锁的机理与应对策略!》此处不再赘述。

剪切自锁仅影响弯曲载荷的完全积分线性单元行为。在受轴向和剪切荷载时,这类单元其实表现很好。

② 二阶单元
从上表还可以看到,不论全积分还是减缩积分的二阶单元表现都很好,自由端的位移结果与理论结果很接近。这是由于二阶单元的边可以弯曲,如下图。但是如果二次单元发生扭曲或弯曲应力有梯度,二次单元也可能发生某种程度的剪切自锁。
 
由以上分析可知,只有我们确信载荷只会在模型中产生很小的弯曲时,才可使用完全积分的线性单元。如果对载荷产生的变形类型有怀疑,则应采用不同类型的单元。

③ 减缩积分的线性单元
从表1可以看到,减缩积分的线性单元随着网格划分越细,越接近理论解答。

但减缩积分的线性单元存在“沙漏”数值问题而过于柔软。由于该类单元只有中部一个积分点,容易形成“零能模式”,产生无意义的结果,即单元变形不产生应变能,形成如下图变形模式:
 

ANSYS中,缩减积分单元软件会采用沙漏控制。在模型中应用的单元越多,对沙漏模式的限制越有效,这说明只要合理地采用细划的网格,线性减缩积分单元就能够给出可接受的结果。

对多数问题,采用线性减缩积分单元的细划网格所产生的误差是在一个可接受的范围之内的。从计算的位移结果可见,当采用这类单元模拟承受弯曲载荷的任何结构时,沿厚度方向上至少采用四个单元才会有比较好的结果。


3.2 弯矩和最大等效应力结果对比
理论计算表明,本例悬臂梁固定约束处的最大弯矩为750N.mm。最大等效应力为71.8MPa。

下面两个表分别是ANSYS关于本例的实体单元不同阶次和积分水平的弯矩和最大等效应力计算结果。

表2 不同单元阶次和积分水平的弯矩结果

表3 不同单元阶次和积分水平的等效应力结果

可以看到,对于弯矩的输出,结果基本一致。只有线性减缩积分单元,输出的结果偏小。这是由于积分点的数量太少引起,当细化网格后,其值越来越接近理论解答。

对于应力的输出,1阶的线性单元结果都比理论值偏小,细化网格后其值也在趋近理论解答,但还是偏差较大。因此,笔者认为对于比较关心应力输出结果时,最好采用二阶单元。

3.3 总结
从本文的分析可以看出,对于一个具体的问题的模拟,如何想得到高精度的结果,正确地选择单元非常重要。从本文的结果看,笔者认为在用ANSYS中做力学有限元分析时,建议:

① 一般分析工作采用采用二阶减缩积分单元,可得到精度较高的位移和应力结果解答,同时也不那耗费计算资源。

② 对于应力集中的部位,采用二阶完全积分单元可以提供应力梯度的最好解答。

除此之外,笔者认为在做有限元分析时,还应尽量减小网格的扭曲,规整的网格可得到更好的结果。


好了,以上就是本期全部内容,如果觉得不错,还请点赞、推荐和转发!对文章有建议,还请在评论区留言交流交流~

来源:薛定谔的Cube
SpaceClaimLS-DYNAWorkbenchAbaqus非线性理论材料控制ANSYS
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2026-02-26
最近编辑:6月前
巴郡撸猫人
硕士 签名征集中
获赞 15粉丝 25文章 73课程 0
点赞
收藏
作者推荐

磨损分析案例前瞻,一文速通ANSYS 接触磨损分析理论!

磨损是指在固体与另一物体接融时,其表面材料逐渐损耗的现象。尽管磨损是一种涉及机械和化学过程的复杂现象,但可以通过将接触表面的各种量与材料损耗关联起来的模型来对其进行近似描述。本公众号在《材料磨损与材料性能关系简析!减小摩擦和磨损有哪些措施?》一文中,简单介绍了材料磨损的一些物理量,并给出了磨损与材料硬度、受载性能的关系。本文介绍一下ANSYS中的磨损模型及磨损理论,作为后续案例基础。--01--ANSYS中的磨损模型1ANSYS中磨损模拟的实现方法在ANSYSAPDL中,在接触分析时,当接触节点被移动到新位置,接触变量(例如,接触压力)会发生变化,当考虑磨损时,需要考虑磨损与接触之间的互相影响。ANSYS通过重新定位接触面上的节点位置来近似计算因磨损导致的材料损耗。节点的新坐标由磨损模型确定。在ANSYS中接触磨损可采用Archard磨损模型(也可采用USERWEAR自定义磨损模型)。CONTA172,CONTA174,andCONTA175接触单元支持磨损模型。要激活接触面的磨损功能,需通过TB、WEAR命令将磨损定义为一个材料模型,并将其分配给接触单元。同时还需要通过TBDATA命令来定义磨损属性。可以通过TBFIELD命令与TBDATA结合使用,来定义温度和/或时间的变化的属性,后续会讲到如何定义。TBTEMP命令也可用于定义温度相关的磨损数据。磨损的实现涉及两个阶段。首先,磨损量通过磨损模型计算。接下来,几何形状将更新以考虑磨损影响。磨损模型根据接触节点处的接触结果计算需要移动多少以及朝哪个方向移动接触节点,以模拟磨损。2Archard磨损模型介绍Archard磨损模型将磨损速率W与接触压力P、滑动速度Vrel和材料硬度H相关联。默认情况下,磨损假设发生在表面的内法线方向,即与接触法线相反的方向。但也可以定义任何所需的磨损方向。Archard磨损模型中,接触节点处的磨损速率W由以下公式给出:其中:K——为磨损系数;H——为材料硬度;P——为接触压力;vrel——为相对滑动速度;m——为压力指数;n——为速度指数。Archard模型由TB、WEAR命令定义,其中TB命令定义时,TBOPT=ARCD。模型所需的材料常数在TBDATA命令中以数据项C1至C4的形式指定。第五个常数C5,可进一步控制Archard模型是如何实现的。常数C6、C7和C8则可用于定义磨损方向的方向余弦。相关材料常数如下:材料常数含义C1磨损系数K,磨损系数通常通过实验确定,因为它依赖于特定的材料配对和工作条件。C2材料硬度H,通过实验测量得到,如布氏硬度、洛氏硬度或维氏硬度等。C3压力指数m,压力指数通常是通过实验数据拟合得到的,它可以反映在不同压力条件下材料的磨损行为,描述磨损体积与接触压力关系的指数。C4速度指数n,通过实验数据来确定,它可以反映在不同滑动速度下材料的磨损特性,描述磨损体积与滑动速度之间的关系。C5可选项:C5=0:默认选项,在磨损计算中使用接触压力。C5=1:在磨损计算中使用节点应力。C5=10:计算接触对接触面积上的磨损增量的平均值。在磨损计算中使用接触压力。C5=11:计算接触对接触面积上的磨损增量的平均值。在磨损计算中使用节点应力。C5=-99:仅用于后处理目的的磨损计算。只计算磨损,不移动接触节点,磨损只是后处理变量,不影响计算方案。C6磨损方向的方位角cosinenx(相对于全局X轴)C7磨损方向的方位角cosineny(相对于全局Y轴)C8磨损方向的方位角cosinenz(相对于全局Z轴)以下是Archard磨损模型的示例命令输入:!材料参数K、H、m、n定义,以下仅为演示K=1E-8!磨损系数H=1000!材料硬度m=1!压力指数n=1!速度指数!定义磨损模型,其中MATID是与接触单元关联的材料ID,workbench中插入命令流时,特定接触面可用CID代替,如下:TB,WEAR,CID,,,ARCD!激活磨损模型,TB命令的TBOPT=ARCDTBDATA,1,K,H,m,n!定义磨损参数TBDATA,6,nx,ny,nz!定义磨损方向关于TB、TBDATA等命令各位置输入参数的含义,可参阅本公众号《ANSYS中关键字、实常数、材料属性和材料模型定义常用命令流!》理解,此处不赘述。3磨损实现的控制选项释义材料的硬度和屈服应力通常密切相关。如果接触单元下的实体单元的材料模型是使用TB命令定义的塑性模型,且Lab=BISO,当前屈服应力的平均值可以用来估计硬度,公式为H=屈服应力/3。要激活此选项,应在TBDATA中输入硬度值(C2)为-99。该选项仅适用于上述提到的材料模型,另外接触单元下的实体单元必须是PLANE182,PLANE183,SOLID185,SOLID186,SOLID187,或SOLID285.在求解的每个子步骤中,软件使用前一子步的屈服应力来估算硬度。因此,从屈服应力更新硬度存在一个滞后一个子步骤的情况,且在求解过程的第一个子步骤中无法获取屈服应力数据。正因如此,若采用此选项,则第一个子步骤中不会发生磨损现象。默认情况下,磨损计算是基于接触压力进行的。如果你在TBDATA命令中对第五个常数(C5)输入值为1,那么磨损计算将基于接触单元下方的实体单元的节点应力,而非接触压力。节点应力用于计算沿接触法向方向的力,而该力值则取代了磨损速率方程中的接触压力。节点应力选项适用于对称接触情况,其中磨损会在两个接触表面均得到模拟。对于非常不同的网格间对称接触而言,接触单元下的实体单元中节点应力的分布往往比接触压力更为平滑。因此,利用节点应力来计算磨损可能会导致磨损更为均匀。如果TBDATA命今中的C5设定为10或11,则磨损增量将被平均分配到接触对的接触区域上,从而使因磨损而导致的总体体积损失与每个节点以不同量磨损(C5=0或1)时因磨损而导致的总体体积损失保持一致。平均磨损增量的计算方式为:其中A表示点处的接触面积,表示接触对中所有接触点的接触面积总和。--02--磨损模型定义此部分磨损定义以ANSYS官方文件中的案例为说明,该模型是,一个半径为30mm的半球形铜环在由钢制成的矩形环上旋转,该矩形环的内半径为50mm,外半径为150mm。半球形环从旋转轴的中心(位于100mm处)开始与平环相触。如下图。具体定义可结合案例模型设置。其中,设置接触面为圆弧线,目标面为与圆弧线接触的直线。1磨损问题的接触设置磨损计算可以采用非对称接触或对称接触方式。使用非对称接触时,模拟仅计算接触面的磨损情况,目标面无磨损。采用对称接触时,模拟将同时计算接触面和目标面的磨损情况。即由于磨损只能在具有接触单元的表面上模拟,非对称接触时仅显示接触面的磨损,对称接触时则同时显示接触面和目标面的磨损。模拟时,注意接触单元有以下设置:①采用增强拉格朗日算法(KEYOPT(2)=0)。执行磨损模拟时,可使用以下接触算法之一:增广拉格朗日或罚函数(KEYOPT(2)=0或1)。使用纯拉格朗日接触算法对磨损进行建模可能会导致收敛问题,不建议使用。②设置接触刚度在每次迭代中更新(KEYOPT(10)=2)。③设置接触检测点的位置为节点,垂直于目标表面(KEYOPT(4)=2)。此处设置原因在于:由于模拟磨损需要重新定位接触节点,因此接触检测点必须位于节点(KEYOPT(4)=1或2),或者可以使用投影法接触(KEYOPT(4)=3)。如下图所示:2定义磨损模型前面提到,磨损模型定义需要通过TB,WEAR命令定义,并赋予接触单元。磨损表面必须被定义为接触单元(而不是目标单元)。在定义Archard磨损模型时,磨损系数K有时可以缩放以简化建模。例如,被磨损件以固定速度旋转时,这种旋转/滑动在接触面上的唯一影响是产生磨损。磨损系数K可以通过建立与滑动速度的关系进行缩放,使得旋转不被显式地建模,但是其影响被包括在磨损的计算中。这将大大减少了模拟时间和工作量。更具体地说,假设磨损率与滑动速度呈线性关系,对于旋转磨损时,滑动速度为2π*转速*旋转半径。将磨损系数K缩放(2π*转速*旋转半径)会使磨损率与滑动速度呈线性关系,而无需对滑动进行明确建模。假设所有点与旋转轴(R)的距离都是恒定的,考虑平均的磨损率,则R取环中心与旋转轴的距离。①单个接触面的磨损(非对称接触)采用非对称接触,仅用于模拟半球形铜环的磨损。在这种情况下,接触单元在铜环上定义,而目标单元在钢环上定义。磨损的材料数据是通过在接触下插入命令流,在TBDATA命令中定义的。设铜环的磨损性能如下:为了在负载达到稳定状态后开始磨损,需分成两个个载荷步进行模拟,定义时TB,wear与TBFIELD,TIME命令需结合使用。在第一个载荷步骤中(TIME=0~1),压力逐渐上升到所需水平,在此载荷步骤中磨损不起作用。为了实现这一目标,WEAR的定义如下:TB,WEAR,CID,,,ARCD!激活接触单元的磨损模型TBFIELD,TIME,0!加载步1开始时的时间TBDATA,1,0,1,1,0,0!设置加载步1开始时磨参!数,C1=0表示无磨损TBFIELD,TIME,1!加载步1结束时的时间TBDATA,1,0,1,1,0,0!设置加载步1结束时磨损!参数,C1=0表示无磨损在第二加载步骤中(TIME=1.01~4),压力保持恒定,发生磨损,磨损定义如下:TBFIELD,TIME,1.01!设置加载步2开始的时间TBDATA,1,kcopper,1,1,0,0!设置加载步2开始时!磨损参数,磨损系数!K=kcopperTBFIELD,TIME,4!设置加载步2结束的时间TBDATA,1,kcopper,1,1,0,0!加载步2期间磨损系!数保持常数K=kcopper②两个接触面的磨损(对称接触)为了模拟两个环上的磨损,需定义对称接触,即两个环上均覆盖一层接触单元,并定义铜环和钢环的磨损特性。铜环的磨损性能如上述非对称所示。钢环的磨损性能如下:不同网格和材料之间的对称接触定义可能会导致接触压力分布不平滑。因此,建议使用接触单元下的实体单元的节点应力来计算磨损增量(TBDATA中定义C5=1)的选项,并在对称示例中使用。铜环的磨损定义如非对称中所示,但C5=1,定义如下:TB,WEAR,CID,,,ARCDTBFIELD,TIME,0TBDATA,1,0,1,1,0,1TBFIELD,TIME,1TBDATA,1,0,1,1,0,1TBFIELD,TIME,1.01TBDATA,1,kcopper,1,1,0,1TBFIELD,TIME,4TBDATA,1,kcopper,1,1,0,1钢环的磨损定义如下:TB,WEAR,TID,,,ARCD!设置为了对称接触,此处可用TID指代目标面。对称接触时,接触面和目标面都覆盖了一层接触单元TBFIELD,TIME,0TBDATA,1,0,1,1,0,1TBFIELD,TIME,1TBDATA,1,0,1,1,0,1TBFIELD,TIME,1.01TBDATA,1,ksteel,1,1,0,1TBFIELD,TIME,4TBDATA,1,ksteel,1,1,0,1③userwear自定义磨损子程序自定义磨损子程序使用了类似于Archard磨损模型的模型。此例中,输入材料数据与Archard模型相似,但增加了角速度。userwear子程序需通过TBDATA定义五个输入属性。C1至C4与Archard磨损定律相同,C5的转速为10000转/秒。示例中的用户磨损子程序使用接触点与旋转轴(R)的距离计算位置相关的滑动速度,并相应地定义磨损增量,从而避免对滑动运动进行显式建模。在这个例子中考虑了不对称接触。铜环的磨损定义如下:V=21e5!将角速度传递给用户进行!计算精确的滑动速度TB,WEAR,CID,,5,USER!激活接触单元自定义磨损TBFIELD,TIME,0!加载步1开始时的时间TBDATA,1,0,1,1,0,V!设置加载步1开始时磨损!参数,C1=0表示无磨损TBFIELD,TIME,1!加载步1结束时的时间TBDATA,1,0,1,1,0,V!设置加载步1结束时磨损!参数,C1=0表示无磨损TBFIELD,TIME,1.01!设置加载步2开始的时间TBDATA,1,kcopper,1,1,0,V!设置加载步2开始时!的磨损参数TBFIELD,TIME,4!设置加载步2结束的时间TBDATA,1,kcopper,1,1,0,V!设置加载步2结束时!的磨损参数3求解过程中改善网格质量磨损建模涉及重新定位接触面节点以模拟材料去除过程。结果,接触单元下方的实体单元的质量会迅速恶化。解决此问题有两种方法,以模拟大量磨损。①非线性网格自适应改进网格的一种方法是使用非线性网格自适应特性。当网格变形时,基于磨损的接触准则会触发非线性网格自适应。磨损量与基础实体单元高度之间的临界比是自定义的。当达到标准时,会触发非线性网格自适应。非对称接触和对称接触都使用非线性网格自适应,需要以下步骤:•创建一个包含正在磨损的接触单元组件;•启用NLADAPTIVE,根据磨损标准触发自适应。假设定义的组件名称是“conwearel”,以下命令将会触发网格自适应:NLADAPTIVE,conwearel,add,contact,wear,0.50!定义磨损标准为单元高度的50%在这种情况下,每当任何接触点的磨损超过接触单元下方实体单元平均高度的50%时,就会发生自适应。每次达到标准时,分析都会停止,通过变形网格来提高网格质量,映射历史相关变量和边界条件,并使用改进的网格重新开始分析,此过程由软件自动完成。关于NLADAPTIVE命令的用法详见《ANSYS中单元类型、材料表字段变量、非线性自适应网格定义命令详解!》一文或帮助文档。②手动重新分区手动重新分区是另一种方法,能够重新对扭曲的网格进行网格划分,并使用改进的网格继续磨损模拟。这种方法需要更多的用户干预。4分析设置问题磨损分析中,存在尺寸的变化,即包括几何非线性,通常需要打开大变形开关。分析设置会使用自动时间增量,磨损过程中接触节点的重新定位可能会导致接触状态的变化。如果磨损增量太大,所有接触单元可能会从闭合状态变为打开状态,导致刚体运动。为了防止这种情况,通常使用非常小的时间增量,这样磨损增量也很小,接触状态的变化也最小化。但这也无疑增加了计算量。总之,使用非常小的子步骤,磨损增量很小,大的磨损增量会突然改变接触状态并导致收敛困难。来源:薛定谔的Cube

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习计划 福利任务
下载APP
联系我们
帮助与反馈