一种基于低质量非结构网格改进的稳定节点积分算法

1.本发明涉及一种基于低质量非结构网格改进的稳定节点积分算法,属于数值算法技术领域。


背景技术：

2.自上个世纪以来，随着电子计算机性能的快速提升，数值计算方法也得到了迅速发展。有限元法(fem)就是以变分原理和分片插值为基础的一种数值计算方法。从简单的静力分析发展到动态、非线性、多物理场耦合等复杂问题的求解，其应用广度与深度都得到了极大的拓展。因其计算效率高，单元划分方式较多，可对形式多样的复杂几何体进行模拟，因此有限元法在实际工程分析中得到了广泛的应用。
3.2007年，刘桂荣等首次将应变光滑技术和有限单元法相结合，提出了光滑有限单元法(sfem)的概念。随后，在总结所提方法特点的基础上，刘桂荣教授提出了g空间理论，由此奠定了一系列方法的理论基础，其中就包含基于节点的光滑有限元法(ns-fem)。
4.ns-fem方法具有很多优良的特性：基于线性三角形或四面体背景网格构建，此种非结构网格可适应复杂的几何结构并可通过商业软件中的自动划分功能很容易地得到；采用线性插值，计算简单、成本低；在ns-fem方法的基础上，后来的研究者又提出了稳定节点积分光滑有限元法(sns-fem)，通过重新构建适宜的积分点来实现在保留节点积分算法特性的同时可以有效地增加系统稳定性并提高计算精度，并且易于实施，还可以有效改善原始方法不稳定、易畸变、精度低的问题。
5.虽然fem方法、ns-fem方法、sns-fem方法已发展的比较全面和成熟，随着研究的不断深入，其在某些问题中的应用仍存在一些不完善之处。比如这三种方法的计算精度都依赖于网格质量和所采用的单元类型。工程分析常用的四边形和六面体单元计算精度较高，但此类网格不能自动生成，在复杂问题分析时需要较高的人工成本。当网格质量降低，不再是所要求的高质量网格，而是不规则的低质量网格(比如较乱的细长、扁平三角形单元网格)时，计算精度会大幅度降低，不稳定性也明显增加。


技术实现要素：

6.为了解决现有技术存在的问题，即fem方法、ns-fem方法、sns方法在低质量网格(比如较乱的细长、扁平三角形单元网格)中计算精度低、不稳定的问题，提供一种基于低质量非结构网格改进的稳定节点积分算法。
7.为了实现本发明的上述目的，本发明提供了一种基于低质量非结构网格改进的稳定节点积分算法，其特殊之处在于包括如下步骤：
8.s1，基于二维低质量非结构网格(非结构网格，通常指的是非四边形和六面体网格，例如三角形网格和四面体网格都属于非结构网格)，输入所需数据，三角形单元的节点编号和节点坐标；
9.s2，构造节点k的光滑域
10.优选的，所述步骤s2中节点k的光滑域构造具体步骤为：
11.在二维情况下，基于低质量(在三角形单元网格中，划分网格质量的标准一般为单元的纵横比，所谓单元纵横比是指单元的高度与宽度之比，最好的值为1，单元形状越不接近等边三角形，三角形越细或者越扁，纵横比也就越大或越小，网格质量越低，本发明所提的低质量网格，三角形单元不但扁或细，排列还较乱，质量更差)的三角形背景网格(比如较乱的细长、扁平三角形单元网格)，依次连接节点k周围三角形单元的面心和边的中点，以得到围绕此节点的多边形作为积分区域ω，积分区域ω被进一步的划分为n
node
个以节点为中心、互不重合的光滑域
12.s3，根据步骤s1中的单元节点编号和节点坐标，对所有节点按照节点编号顺序循环，判断节点是在内部还是在边界上，节点若在内部，将步骤s2中节点k的多边形光滑域近似为等面积的椭圆形；节点若在边界上，将步骤s2中节点k的多边形光滑域近似为等面积的圆形；
13.s4，根据步骤s3中的已知条件，求出椭圆形光滑域的长短半轴lc、ld，或者是圆形光滑域的半径lr，以及与光滑域内局部坐标系坐标轴的夹角α；
14.优选的，所述步骤s4中节点在内部和边界上时的长度和夹角具体为：
15.基于低质量网格，将节点周围的多边形光滑域近似为等面积的椭圆形，作为近似光滑域当节点为单元内部点时，以多边形上离节点(xk,yk)最远距离的一个点(x0,y0)作为近似椭圆长轴的一个端点，并以此最远距离为椭圆长半轴长度lc，根据椭圆面积与多边形面积相等的条件即可算出椭圆短半轴长度ld，其中lc、ld的计算公式(1)、(2)为：
[0016][0017][0018]
式中，ak为多边形面积，然后将进一步划分为4个子光滑域，选择椭圆长轴和短轴的四个端点为积分点并取每个子光滑域内局部坐标系坐标轴与椭圆长、短轴的夹角α，夹角α的计算公式(3)为：
[0019][0020]
当节点为边界点时，将节点周围的多边形光滑域近似为等面积的圆形，此时椭圆形近似光滑域的长短半轴变为同一长度lr，即圆形近似光滑域的半径，且与光滑域内局部坐标系坐标轴的夹角α＝0；
[0021]
s5，求出四个积分点处的应变
[0022]
所述步骤s5中四个积分点的应变具体为：
[0023]
有了积分点和积分域之后，应变也可以通过计算得到，假定应变在内是连续且一阶可导的，其泰勒展开为公式(4)：
[0024][0025]
式中，εk表示k节点处的应变，四个子域积分点处的应变表示k节点处的应变，四个子域积分点处的应变分别为公式(5)-(12):
[0026][0027][0028][0029][0030]
l
d1
＝ld·
sinα
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
[0031]
l
d2
＝ld·
cosα
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0032]
l
c1
＝lc·
sinα
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(11)
[0033]
l
c2
＝lc·
cosα
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(12)
[0034]
s6，由步骤s5中的应变得出应变能；
[0035]
所述步骤s6中的应变能为四个积分子域上应变能，应变能为公式(13)：
[0036][0037]
其中，d为弹性模量的四阶张量，将代入上式得：
[0038][0039]
式中，
[0040][0041]
s7，根据步骤s6中的应变能求出每个节点积分域的稳定刚度矩阵
[0042]
所述步骤s7中的稳定刚度矩阵根据变分原理，刚度矩阵可按下式公式(14)计算：
[0043][0044]
式中，k
sc
为总体的稳定刚度矩阵，为每个节点积分域的稳定刚度矩阵，的计算公式(15)为:
[0045][0046]
式中bk，δb
xy1
，δb
xy2
为应变转换矩阵，均采用梯度光滑技术将区域积分转化为边界积分计算，其中：
[0047][0048][0049][0050]
式中，表示形函数一阶导数，n
i,jj
(jj＝xx,xy,yx,yy)表示形函数二阶导数，,其中表示积分域内的节点数目；
[0051]
s8，根据总的刚度矩阵算出节点位移和节点应力结果。
[0052]
本发明通过改进的方法求出刚度矩阵之后，然后再根据总的刚度矩阵得到更加精确的节点位移和节点应力结果。与现有技术相比，本发明的有益效果在于：
[0053]
提出了一种在低质量非结构网格中具有高精度、高收敛的稳定节点积分算法，在保证计算效率的前提下，极大改善了算法在低质量网格的计算精度。
附图说明
[0054]
图1是基于低质量非结构网格改进的稳定节点积分算法的流程图；
[0055]
图2是本发明一种优选实施方式中采用的二维低质量网格问题中基于节点k的光滑域示意图；
[0056]
图3是本发明一种优选实施方式中椭圆近似积分域内部具体构造图；
[0057]
图4是本发明一种优选实施方式中二维静态问题中的结构示意图；
[0058]
图5是本发明一种优选实施方式中悬臂梁的低质量三角形网格划分示意图；
[0059]
图6a是本发明一种优选实施方式中的不同算法在低质量网格下的悬臂梁中心线位移结果的比较示意图；
[0060]
图6b是本发明一种优选实施方式中的不同算法在低质量网格下的悬臂梁右端中心点位移收敛曲线的比较示意图；
[0061]
图6c是本发明一种优选实施方式中的不同算法在低质量网格下的悬臂梁系统应变能收敛曲线的比较示意图。
[0062]
图7是本发明一种优选实施方式中二维自由振动问题中的简化车门的结构示意图；
[0063]
图8是本发明一种优选实施方式中简化车门的低质量三角形网格划分示意图；
[0064]
图9是本发明一种优选实施方式中的不同算法在低质量网格下的简化车门前24阶自由振动固有频率相对误差收敛结果的比较示意图。
具体实施方式
[0065]
为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明的技术方案进行清楚、完整地描述。本技术所描述的实施例仅仅是本发明一部分的实施例，而不是全部实施例。基于本发明精神，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其它实施例，都属于本发明的保护范围。
[0066]
实施例1
[0067]
本发明提供了一种基于低质量非结构网格改进的稳定节点积分算法，如附图1-3所示，下面结合附图1-3对本发明的改进方法进行详细说明。
[0068]
本发明的基于低质量非结构网格改进的稳定节点积分算法，包括如下步骤：
[0069]
s1，基于二维低质量非结构网格(比如较乱的细长或扁平三角形网格)，输入所需数据，三角形单元的节点编号和节点坐标；
[0070]
s2，构造节点k的光滑域参考附图2显示了二维问题中基于节点的光滑域构造方法，光滑域的边界用表示，在二维情况下，基于低质量的三角形背景网格，依次连接节点k周围三角形单元的面心和边的中点,以得到围绕此节点的多边形作为积分区域ω，积分区域ω被进一步的划分为n
node
个以节点为中心、互不重合的光滑域
[0071]
s3，根据步骤s1中的单元节点编号和坐标，对所有节点按照节点编号顺序循环，判断节点是在内部还是在边界上，节点若在内部，将步骤s2中节点k的多边形光滑域近似为等面积的椭圆形；节点若在边界上，将步骤s2中节点k的多边形光滑域近似为等面积的圆形。
[0072]
s4，根据步骤s3中的已知条件，求出椭圆形光滑域的长短半轴lc、ld，或者是圆形光滑域的半径lr，以及与光滑域内局部坐标系坐标轴的夹角α。基于低质量网格，当节点为单元内部点时，将节点周围的多边形光滑域近似为等面积的椭圆形，作为近似光滑域如图3所示。以多边形上离节点(xk,yk)最远距离的一个点(x0,y0)作为近似椭圆长轴的
一个端点，并以此最远距离为椭圆长半轴长度lc，根据椭圆面积与多边形面积相等的条件即可算出椭圆短半轴长度ld的计算公式为：
[0073][0074][0075]
式中，ak为多边形面积，然后将进一步划分为4个子光滑域，选择椭圆长轴和短轴的四个端点为积分点并取每个子光滑域内局部坐标系坐标轴与椭圆长、短轴的夹角α，如图3所示。夹角α的计算公式为：
[0076][0077]
当节点为边界点时，将节点周围的多边形光滑域近似为等面积的圆形，此时椭圆形近似光滑域的长短半轴变为同一长度lr，即圆形近似光滑域的半径，且与光滑域内局部坐标系坐标轴的夹角α＝0；
[0078]
s5，求出四个积分点处的应变有了积分点和积分域之后，应变也可以通过计算得到。假定应变在内是连续且一阶可导的，其泰勒展开为：
[0079][0080]
式中，εk表示k节点处的应变。四个子域积分点处的应变表示k节点处的应变。四个子域积分点处的应变分别为:
[0081][0082][0083][0084][0085]
l
d1
＝ld·
sinα
[0086]
l
d2
＝ld·
cosα
[0087]
l
c1
＝lc·
sinα
[0088]
l
c2
＝lc·
cosα
[0089]
s6，由步骤s5中的应变得出应变能：
[0090][0091]
其中，d为弹性模量的四阶张量，将代入上式可得：
[0092][0093]
式中，
[0094][0095]
s7，根据步骤s6中的应变能求出每个节点积分域的稳定刚度矩阵根据变分原理，刚度矩阵可按下式计算：
[0096][0097]
式中，k
sc
为总体的稳定刚度矩阵，为每个节点积分域的稳定刚度矩阵，的计算公式为:
[0098][0099]
式中bk，δb
xy1
，δb
xy2
为应变转换矩阵，均采用梯度光滑技术将区域积分转化为边界积分计算，其中:
[0100][0101][0102][0103]
式中，表示形函数一阶导数,n
i,jj
(jj＝xx,xy,yx,
sns-fem)与有限元法(fem)、基于节点的光滑有限元法(ns-fem)、稳定节点积分光滑有限元法(sns-fem)的前24阶自由振动固有频率相对误差结果一起绘制在附图9中进行对比，可以发现本发明的算法(m-sns-fem)的总体相对误差结果更加接近于0，且所有误差均在正负10％以内(正负10％为相对误差允许的最大误差范围，若超过这个范围值则说明精度很低)，说明本发明的算法误差最小，比其它方法更加稳定和精确。
[0119]
本发明改进的稳定节点积分算法在不增加计算成本的前提下，降低了对网格质量的要求，极大改善了有限元的计算精度，提高了算法的适用性。
[0120]
以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等同替换，其均应涵盖在本发明的权利要求保护范围之内。