有限元的误差主要来自两个方面:一是模型误差,一是计算误差。
模型误差:是指将实际工程问题抽象为适合计算机求解的有限元模型时所产生的误差,即有限元模型和实际问题之间的差异。它包括有限元离散处理所固有的原理性误差,也包括几何模型处理、实际工况转化为模型边界条件时所带来的偶然性误差。
几何离散误差是指离散后的几何体与原有几何形状上的差异。对于由直线或者平面边界构成的规则结构这类误差较小。对于具有复杂曲线或者曲面边界构成的结构离散后会产生较大的形状误差。
(离散误差 : 我们都知道目前的有限元分析方法是基于位移法的,而位移的表述我们自然要用到形函数。(常规的就是三角形单元,四边形单元的形函数)
物理离散误差是插值函数和真实函数之间的差异,其大小与单元尺寸和插值多项式的阶次有关。一般单元尺寸减小也就是网格划分越密,插值函数的阶次增加,将使有限元的解收敛于精确解。
计算误差:是指采用数值方法对有限元模型进行计算所产生的误差,误差的性质是舍入误差和截断误差。
那么,显而易见,有限元计算的精度与模型误差和计算误差存在着极大的关系,其中几何模型的网格划分层次和区域,网格数量和质量,以及单元类型的选取起着至关重要的作用。
选择单元首先需要明确单元的类型,结构有限元分析中主要有以下一些单元类型:平面应力单元、平面应变单元、轴对称实体单元、空间实体单元、板单元、壳单元、 轴对称壳单元、杆单元、梁单元、弹簧单元、间隙单元、质量单元、摩擦单元、刚体单元和约束单元等。根据不同的分类方法,上述单元可以分成以下不同的形式。
按照维度进行单元分类 :
根据单元的维数特征单元可以分为一维单元、二维单元和三维单元。一维单元的网格为一条直线或者曲线。直线表示由两个节点确定的线性单元。曲线代表由两个以上的节点确定的高次单元,或者由具有确定形状的线性单元。杆单元、梁单元.
二维单元的网格是一个平面或者曲面它没有厚度方向的尺寸。这类单元包括平面单元、轴对称实体单元、板单元、壳单元和复合材料壳单元等。二维单元的形状通常具有三角形和四边形两种。在使用自动网格剖分时,这类单元要求的几何形状是表面模型或者实体模型的边界面,或者抽取中面。采用薄壳单元通常具有相当好的计算率。
三维单元的网格具有空间三个方向的尺寸,其形状具有四面体、五面体和六面体。这类单元包括空间实体单元和厚壳单元。在自动网格划分时,它要求的是几何模型是实体模型(厚壳单元是曲面也可以)。
按照插值函数进行单元分类:
根据单元插值函数多项式的最高阶数多少,单元可以分为线性单元、二次单元、三次单元和更高次的单元。
线性单元具有线性形式的插值函数,其网格通常只具有角节点而无边节点,网格边界为直线或者平面。这类单元的优点是节点数量少,在精度要求不高或者结果数据梯度不太大的情况下,采用线性单元可以得到较小的模型规模。注意:但是由于单元位移函数是线性的,单元内的位移呈线性变化,而应力是常数。因此会造成单元间的应力不连续,单元边界上存在着应力突变。
二次单元的插值函数是二次多项式其网格不仅在每个顶点处有角节点,而且在棱边上还存在一个边节点,因此网格边界可以是二次曲线或者曲面。这类单元的优点是几何和物理离散精度较高,单元内的位移呈二次变化,应力呈线性变化,因此单元边界上的应力是连续的。但是在单元数量相同的条件下二次单元的节点数比线性单元的节点数多,模型的规模较大。
三次单元的插值函数是三次多项式,其网格的每条边上存在两个节点。有些三次单元还具有内部节点。这类单元的离散精度更高,但是由于单元节点数较多,网格划分较为困难,模型规模很大,一般用于具有特殊精度要求的场合。
在有限元分析中,一般认为当节点数目或单元插值位移趋于无穷大时,即当单元尺寸趋近于零时,最后的解答如果能够无限地逼近准确解,那么这样的位移函数(形状函数)是逼近于真解的,这就称为收敛。
右图表示出几种可能的收敛情况。
很明显,其中曲线1和曲线2都是收敛的,但曲线1比曲线2收敛更快。曲线3虽然趋向于某一确定值,但该值不是问题的准确解,故也不能算作收敛。曲线4虽然收敛,但不是单调收敛,所以也不能算是收敛。至于曲线5,则是发散,完全不符合要求。
通过以上描述,我们注意到几个关键词:节点数目,单元插值位移 ,形状函数,收敛更快,趋近于某一确定值…
网格数量的多少将影响计算结果的精度和计算规模的大小。一般来讲,网格数量增加,计算精度会有所提高,但同时计算规模也会增加,所以在确定网格数量时应权衡两个因数综合考虑。
网格较少时增加网格数量可以使计算精度明显提高,而计算时间不会有大的增加。当网格数量增加到一定程度后,再继续增加网格时精度提高甚微,而计算时间却有大幅度增加。所以应注意增加网格的经济性。实际应用时可以比较两种网格划分的计算结果,如果两次计算结果相差较大,可以继续增加网格,相反则停止计算。
在决定网格数量时应考虑分析数据的类型。在静力分析时,如果仅仅是计算结构的变形,网格数量可以少一些。如果需要计算应力,则在精度要求相同的情况下应取相对较多的网格。同样在响应计算中,计算应力响应所取的网格数应比计算位移响应多。在计算结构固有动力特性时,若仅仅是计算少数低阶模态,可以选择较少的网格,如果计算的模态阶次较高,则应选择较多的网格。在热分析中,结构内部的温度梯度不大,不需要大量的内部单元,这时可划分较少的网格。
网格疏密是指在结构不同部位采用大小不同的网格,这是为了适应计算数据的分布特点。在计算数据变化梯度较大的部位(如应力集中处),为了较好地反映数据变化规律,需要采用比较密集的网格。而在计算数据变化梯度较小的部位,为减小模型规模,则应划分相对稀疏的网格。这样,整个结构便表现出疏密不同的网格划分形式。采用疏密不同的网格划分,既可以保持相当的计算精度,又可使网格数量减小。因此,网格数量应增加到结构的关键部位,在次要部位增加网格是不必要的,也是不经济的。
划分疏密不同的网格主要用于应力分析(包括静应力和动应力),而计算固有特性时则趋于采用较均匀的网格形式。这是因为固有频率和振型主要取决于结构质量分布和刚度分布,不存在类似应力集中的现象,采用均匀网格可使结构刚度矩阵和质量矩阵的元素不致相差太大,可减小数值计算误差。同样,在结构温度场计算中也趋于采用均匀网格。
许多单元都具有线性、二次和三次等形式,其中二次和三次形式的单元称为高阶单元。选用高阶单元可提高计算精度,因为高阶单元的曲线或曲面边界能够更好地逼近结构的曲线和曲面边界,且高次插值函数可更高精度地逼近复杂场函数,所以当结构形状不规则、应力分布或变形很复杂时可以选用高阶单元。
但高阶单元的节点数较多,在网格数量相同的情况下由高阶单元组成的模型规模要大得多,因此在使用时应权衡考虑计算精度和时间。在网格数量较少时,两种单元的计算精度相差很大,这时采用低阶单元是不合适的。当网格数量较多时,两种单元的精度相差并不很大,这时采用高阶单元并不经济。
(例如在离散细节时,由于细节尺寸限制,要求细节附近的网格划分很密,这时采用线性单元更合适。 增加网格数量和单元阶次都可以提高计算精度。因此在精度一定的情况下,用高阶单元离散结构时应选择适当的网格数量,太多的网格并不能明显提高计算精度,反而会使计算时间大大增加。)
为了兼顾计算精度和计算量,同一结构可以采用不同阶次的单元,即精度要求高的重要部位用高阶单元,精度要求低的次要部位用低阶单元。不同阶次单元之间或采用特殊的过渡单元连接,或采用多点约束等式连接。
当然,对于特定的分析也应选用特定阶次的单元,(例如:对于变形体之间的通用接触分析,最佳使用一阶四边形或者六面体单元,应尽量避免使用二次四边形和六面体单元,至于单元选取,这个我们后边的课程会讲到。)
工程中我们经常会推荐大家使用四边形或者六面体单元。这里我们将对三节点三角形单元与四节点四边形单元的计算精度做以比较,以说明选取这类单元的必要性。
首先我们从单元的形函数来看:
三角形单元的位移场为线性关系式,平面位移模式如下:u(x,y)=a0+a1x+a2y,v(x,y)=b0+b1x+b2y)
求解组成矩阵为常系数矩阵,单元内任意一点的应变应力都是常数。因此,三节点三角形单元称为常应力(应变)CST单元(constant strain triangle).在实际使用过程中,对于应变(应力)梯度较大的区域,单元划分应适当紧密,否则将不能反映应变或者应力的真实变化情况,从而导致较大的误差。
四边形单元的位移场为双线性位移模式,平面位移模式如下:
u(x,y)=a0+a1x+a2y+a3xy,v(x,y)=b0+b1x+b2y+b3xy
在单元的边界x=+(-)a和y=+(-)b上,位移是按照线性变化的,且相邻单元公共节点上具有共同的节点位移值,可保证两个相邻单元在其公共边界上的位移是连续的,这种单元的位移模式是完备和协调的,他的应力应变呈现一次性变化,因而比三节点常应变单元精度高。
位移场和应力场:
同一工况下的两种单元精度对比:Mises应力结果(外插平均后)
同一工况下的两种单元对比:位移结果对比(外插平均后)
3节点三角形单元与四节点矩形单元的计算精度比较:
网格质量是指网格几何形状的合理性。质量好坏将影响计算精度。质量太差的网格甚至会中止计算。直观上看,网格各边或各个内角相差不大、网格面不过分扭曲、边节点位于边界等份点附近的网格质量较好。网格质量可用细长比、锥度比、内角、翘曲量、拉伸值、边节点位置偏差等指标度量。
划分网格时一般要求网格质量能达到某些指标要求。在重点研究的结构关键部位,应保证划分高质量网格,即使是个别质量很差的网格也会引起很大的局部误差。而在结构次要部位,网格质量可适当降低。当模型中存在质量很差的网格(称为畸形网格)时,计算过程将无法进行。
结构中的一些特殊界面和特殊点应分为网格边界或节点以便定义材料特性、物理特性、载荷和位移约束条件。即应使网格形式满足边界条件特点,而不应让边界条件来适应网格。
常见的特殊界面和特殊点有材料分界面、几何尺寸突变面、分布载荷分界线(点)、集中载荷作用点和位移约束作用点等。
当结构形状对称时,其网格也应划分对称网格,以使模型表现出相应的对称特性(如集中质矩阵对称)。不对称布局会引起一定误差
对于具有旋转对称,轴对称的有限元模型,都应按照对称方式来做以求达到更高的精度,这些内容在前处理软件中很容易做到。
节点个数 (插值)
节点的个数决定了单元的插值方式。
ABAQUS包含一阶和二阶插值方式的单元。
所谓“完全积分”是指当单元具有规则形状时,所用的Gauss积分点的数目足以对单元刚度矩阵中的多项式进行精确积分。对六面体和四边形单元而言,所谓“规则形状”是指单元的边是直线并且边与边相交成直角,在任何边中的节点都位于边的中点上。完全积分的线性单元在每一个方向上采用两个积分点。因此,三维单元C3D8在单元中采用了222个积分点。完全积分的二次单元(仅存在于ABAQUS/Standard)在每一个方向上采用3个积分点。对于二维四边形单元,完全积分的积分点位置如图所示。
对于线性完全积分单元,在厚度方向的单元数目并不影响计算结果。误差是由于剪力自锁(shear locking)引起的,这是存在于所有完全积分、一阶实体单元中的问题。
剪力自锁引起单元在弯曲时过于刚硬,对此解释如下。考虑受纯弯曲结构中的一小块材料,如图4-4所示,材料产生弯曲,变形前平行于水平轴的直线成为常曲率的曲线,而沿厚度方向的直线仍保持为直线,水平线与竖直线之间的夹角保持为90度 。
线性单元的边不能弯曲;所以,如果应用单一单元来模拟这一小块材料,其变形后的形状如图所示。
受弯矩M作用下完全积分、线性单元的变形
为清楚起见,画出了通过积分点的虚线。显然,上部虚线的长度增加,说明1方向的应力是拉伸的。类似地,下部虚线的长度缩短,说明是压缩的。竖直方向虚线的长度没有改变(假设位移是很小的);所有这些都与受纯弯曲的小块材料应力的预期状态是一致的。但是,在每一个积分点处,竖直线与水平线之间夹角开始时为90度,变形后却改变了,说明这些点上的剪应力不为零。显然,这是不正确的:在纯弯曲时,这一小块材料中的剪应力应该为零。
产生这种伪剪应力的原因是因为单元的边不能弯曲,它的出现意味着应变能正在产生剪切变形,而不是产生所希望的弯曲变形,因此总的挠度变小:即单元是过于的刚硬。
剪力自锁仅影响受弯曲载荷的完全积分的线性单元的行为。在受轴向或剪切载荷时,这些单元的功能表现很好。而二次单元的边界可以弯曲,故它没有剪力自锁的问题。二次单元预测的自由端位移接近于理论解答。
受弯矩M作用下完全积分、二次单元的变形
只有当确信载荷只会在模型中产生很小的弯曲时,才可以采用完全积分的线性单元。
如果对载荷产生的变形类型有所怀疑,则应采用不同类型的单元。
在复杂应力状态下,完全积分的二次单元也有可能发生自锁;因此,如果在模型中应用这类单元,应细心地检查计算结果。
然而,对于模拟局部应力集中的区域,应用这类单元是非常有用的!
减缩积分
只有四边形和六面体单元才能采用减缩积分方法;
而所有的楔形体、四面体和三角形实体单元采用完全积分,尽管它们与减缩积分的六面体或四边形单元可以在同一网格中使用。
减缩积分单元比完全积分单元在每个方向少用一个积分点。减缩积分的线性单元只在单元的中心有一个积分点。(实际上,在ABAQUS中这些一阶单元采用了更精确的均匀应变公式,即计算了单元应变分量的平均值。对于所讨论的这种区别并不重要。)对于减缩积分的四边形单元,积分点的位置如图所示。
线性的减缩积分单元由于存在着来自本身的所谓沙漏(hourglassing)数值问题而过于柔软。为了说明这个问题,再次考虑用单一减缩单元模拟受纯弯曲载荷的一小块材料 :
受弯矩M的减缩积分线性单元的变形
单元中虚线的长度没有改变,它们之间的夹角也没有改变,这意味着在单元单个积分点上的所有应力分量均为零。由于单元变形没有产生应变能,这种变形的弯曲模式是一个零能量模式。由于单元在此模式下没有刚度,所以单元不能抵抗这种形式的变形。在粗划网格中,这种零能量模式会通过网格扩展,从而产生无意义的结果。
ABAQUS在一阶减缩积分单元中引入了一个小量的人工“沙漏刚度”以限制沙漏模式的扩展。在模型中应用的单元越多,这种刚度对沙漏模式的限制越有效,这说明只要合理地采用细划的网格,线性减缩积分单元可以给出可接受的结果。对多数问题而言,采用线性减缩积分单元的细划网格所产生的误差(见表4-2)是在一个可接受的范围之内。结果建议当采用这类单元模拟承受弯曲载荷的任何结构时,沿厚度方向上至少应采用四个单元。当沿梁的厚度方向采用单一线性减缩积分单元时,所有的积分点都位于中性轴上,该模型是不能抵抗弯曲载荷的。(这种情况在表4-2中用*标出。)
线性减缩积分单元能够很好地承受扭曲变形;因此,在任何扭曲变形很大的模拟中可以采用网格细划的这类单元。
在ABAQUS/Standard中,二次减缩积分单元也有沙漏模式。然而,在正常的网格中这种模式几乎不能扩展,并且在网格足够加密时不会产生什么问题。因此,除了包含大应变的大位移模拟和某些类型的接触分析之外,这些单元一般是最普遍的应力/位移模拟的最佳选择。