基于平均值等几何法的平面静力学分析

孙振山，毛虎平，贾冰潮，张稣宁

(中北大学 能源与动力工程学院，山西 太原 030051)

有限元法是求解静力学问题常用的方法，在平面结构分析中也不例外. 平面静力学分析以力学原理和数值分析算法为基础，通过构造整体刚度矩阵建立本构方程，在指定的设计要求和约束条件下，计算其位移、应力等，进而为实际的生产加工提供参考.

有限元分析法是机械、土木、汽车、船舶及航空航天等领域中平面静力学分析普遍采用的基础手段. 对于结构复杂且需要精确表示的平面，目前依赖于传统的建模方法，利用CAD与CAE建模，然后将模型导入有限元分析软件进行计算. 随着设计模型复杂程度的增加，对有限元网格划分的质量、计算速度、仿真精度有了更高的要求. 在有限元分析过程中网格划分占据约80%的时间，而且有限元分析所建立的几何模型是近似的而非精确的表示，所引入的几何误差对一些关键部件的影响较大. 通常，工程师会通过不断细化网格来达到所需要的计算精度，但这样不仅会产生大量的网格单元，也会降低计算速度，所需要的计算成本也会随之增加. 有限元分析和设计的过程需要计算机不断地修改和重新划分网格，此过程需要与CAD模型频繁交互，过程繁琐. 此外，在有限元分析过程中以C0连续性为代表的常规单元无法满足高精度高性能的分析要求.

研究者希望超越现有网格划分的分析方法，直接基于精确的几何模型进行建模、分析和优化，将所需要的几何模型通过控制点精确地表示出来，但此方法相较于成熟的有限元分析方法而言也存在一些缺点. 2005年，Hughes等[1]提出了基于NURBS(Non-Uniform Rational B-Splines)基函数的等几何分析(Isogeometric Analysis, IGA)数值计算方法，该方法将所需要的几何模型通过NURBS基函数建立精确的参数化模型，通过对节点、节点矢量、控制点及控制点权重的计算，选择合适的分段多项式基函数，构造出分段有理多项式曲线以及非有理的张量积曲面.

等几何分析方法以样条(如B样条、T样条等)函数为基础，将CAD几何模型与结构分析紧密结合，从而实现精确的几何分析[2-3]. 将等几何分析应用于平面静力学分析时，直接建立起了IGA与平面边界的联系，这有利于模型参数化[4-5]. 大部分的CAD和CAE软件在建模时采用布尔和运算，这将会产生大量的剪切曲面，使得设计难度增加，而等几何分析直接基于物理模型建模，省去了大量的剪切操作. 与传统的有限元分析相比，等几何分析在优化过程中有如下优点：(1) 通过精确的NURBS基函数对模型进行建模分析，从源头上消除了有限元分析所引入的几何误差，使得计算精度更高；
(2) 避免了在有限元建模过程中与计算机进行繁琐交互的问题；
(3) 等几何分析将计算所得的控制点作为设计变量，在保证分析精度要求的同时具有更好的连续性，从而有效克服了有限元分析方法中C0连续性所引起的计算困难[6-8]的问题. 目前，国内外学者对等几何分析也开展了大量相关研究，如袁沛等[9]采用了一种基于等几何分析的非结构化T样条构建复杂模型的方法；
贾悦等[10]结合传统等几何分析配置法和传统等几何分析伽辽金法，提出了一种基于PHT-样条函数的加强等几何分析配置方法；
张中林等[11]引入了基于推广B样条的等几何分析方法对二维弹性问题进行边界元分析；
胡丹丹等[12]针对等几何分析中复杂物理区域上的偏微分方程求解问题，提出了多片参数域上双三次样条投影映射的方法. 在等几何法应用方面，李明超等[13]基于等几何分析开展了水利结构数值仿真计算新方法的研究；
柏硌等[14]为解决零亏格边界网格模型的T样条实体重建问题，提出了一种基于八叉树细分和渐进迭代最小二乘拟合算法的T样条实体构建算法；
陈龙等[15]将等几何分析应用于单齿啮合的齿轮接触问题上，此法在解决接触仿真分析方面具有优良的性能.

综上所述，等几何分析由于精确的参数化建模，使得其分析更加方便，求解效率更高，所需要的计算成本相对更低，因此已被广泛应用于许多领域. 但与成熟的有限元技术相比，等几何法起步较晚，目前还没有一款商用软件将等几何分析成功应用于多种工程实例. 节点向量是等几何分析的关键要素之一，节点向量的等距分布可能会产生奇异方程组，降低矩阵计算效率. 本文借助matlab平台，将平均值法融入到节点向量的计算中，并将平均值等几何分析方法应用于经典力学问题，建立了平面静力学系统化程序. 采用平均值方法对节点向量进行计算，以NURBS基函数对平面结构进行参数化建模，以二维平面静力学问题为例展现此算法的可行性.

1.1 NURBS基函数

非均匀有理B样条(NURBS，Non-Uniform Rational B-Splines)以B样条方法为基础，具有非均匀性、有理性和B样条三个特征，即控制点的影响范围可以改变、有理多项式定义几何结构和控制多边形或控制多面体来构造曲线曲面.k次有理基函数

(1)

式中：ωi为对应控制点的权值；
Ni,k(ξ)为k次B样条基函数，由de Boor-Cox[16-18]提出的递推公式定义为

(2)

为了表示面，二维NURBS基函数可表示为

(3)

NURBS曲面是一个双线性分段有理函数，表示为

(4)

式中：Pi,j为控制点.

1.2 节点向量计算方法

建立一个系数矩阵为(n+1)×(n+1)的线性方程组

(5)

式中：{Qk}为一组数据点.

(6)

(7)

u0=…=up=0,um-p=…=um=1,

(8)

式中：n为控制点的最大下标(控制点的个数为n+1)；
p为曲线的次数；
m+1为节点个数(控制点、节点从0开始计数).

1.3 等几何分析

图1 所示为NURBS曲面等几何分析示意图，可以看出，等几何分析首先通过控制点建立物理空间，控制点部分在曲面上，部分不在曲面上，然后采用等参元思想将物理空间转化到[-1，1]的父空间，再将父空间转化到[0，1]参数空间，得到节点向量，最后建立等几何分析模型.

图1 NURBS曲面等几何分析示意图[2]Fig.1 Schematic diagram of geometric analysis of NURBS surface

采用NURBS基函数离散时，涉及到参数变换问题. NURBS基函数定义在标准[0,1]空间，而物理空间不在此空间，因此需要进行转换. 数值积分一般采用Gauss积分，而其积分空间在标准[-1,1]空间，需要从[-1,1]空间转换到[0,1]空间.

(9)

对应的雅可比公式为

(10)

从物理空间[a,b]到标准[0,1]空间转化的雅可比公式为

(11)

在优化中，用等几何分析法代替有限元法对物理场进行数值计算. 采用NURBS基函数可以直接表示几何模型并将其作为结构分析的形函数. 参数坐标(ξ,η)中的物理量可以通过控制点的值获得

(12)

对于一个线弹性问题，控制方程离散后获得

Ku=F,

(13)

式中：K表示全局刚度矩阵;u表示位移向量;F表示与控制点相关的外力.等几何分析中，同样通过单元刚度矩阵Ke来构造全局刚度矩阵K，Ke表示为

(14)

(15)

而

(16)

式中：Ni为NURBS单元的基函数;J1为

(17)

雅可比矩阵J2表示为

(18)

对于完全弹性的各向同性材料，弹性矩阵D表示为

(19)

式中：E为弹性模量；
μ为泊松比.

将单元刚度矩阵Ke逐个按索引编号组装到整体矩阵中，组装过程与传统的有限元方法类似，得到关于所有控制点变量的整体方程

Ku=F.

(20)

整体刚度矩阵

(21)

整体等效控制点载荷向量

(22)

等几何分析与有限元分析类似，采取一种类等参元的思想. 不同之处是等几何分析一般采用NURBS基函数，而有限元法采用拉格朗日基函数. NURBS基函数具有局部支撑性、规范性、可微性和非负性等性质. 其中，可微性使等几何分析建立的曲线曲面更加光滑，连续性更好，同时它又可以出现尖点. NURBS基函数在节点区间内部存在任意阶导数，在节点处，它是p-k次连续可微的，k为该节点的重复度. 等几何分析在建模过程中直接基于精确的几何模型建模，可以避免有限元分析在划分网格时可能引入的误差. 图2 为二维平面静力学等几何分析的基本流程图.

图2 等几何分析平面静力学分析流程图Fig.2 Flow chart of plane statics analysis for isogeometric analysis

本文采用matlab编程实现二维平面问题的等几何分析. 平面问题的等几何分析属于线性静态分析，首先定义平面的材料属性，然后固定该平面的一端使其x轴位移为0，对另一端施加集中力载荷. 之后将平面模型导入，将平面模型从物理域向几何域转化，并参数化几何域模型. 根据算法计算平面模型的节点向量以及基函数，再根据要求对其模型进行节点插入、网格细分. 引入等参元的思想，计算每个等几何单元的刚度矩阵，再将它们整合到整体刚度矩阵，把边界条件施加到整体刚度矩阵中，最终求得平面的位移、应力和应变.

算法具体实现过程如下：

1) 采用式(8)计算节点向量；

2) 定义基本变量，包括控制点、基函数次数、节点向量和材料属性,调用弹性矩阵算法；

3) 进行单元细化，调用单元细化算法(hRefinement2d)和生成等几何二维网格算法(generateIGA2DMesh)；

4) 提取控制点数和单元数量，为节点编号；

5) 通过调用plotPlate算法绘制几何模型；

6) 找寻相应的边界节点编号和受力单元节点编号对其进行初始化设置；

7) Gauss积分；

8) 编写循环单元刚度矩阵，并将之组装到整体刚度矩阵；

9) 对受力点施加节点力F；

10) 进行边界条件处理，简化刚度矩阵；

11) 通过u=KF求解位移；

12) 后处理，绘制位移云图.

为了验证改进后的等几何法的可行性，选择一个普通的平面静力学模型进行计算.

3.1 悬臂梁问题

悬臂梁结构如图3 所示，悬臂梁长L=30 m，梁的高度h=6 m，单位厚度. 平面左端固定，右端施加集中力载荷F=-10 000 N, 方向垂直向下. 材料为线弹性，弹性模量E=210 GPa，泊松比μ=0.3.

图3 悬臂梁结构Fig.3 Cantilever beam structure

3.2 悬臂梁解析解

图4 悬臂梁模型Fig.4 Cantilever beam model

此问题的位移解析解和应力解析解[19]为

u(x,y)=

(23)

v(x,y)=

(24)

(25)

由式(23)～式(25)计算得，此问题x方向的最大位移解析解 (ux)max=3.571×10-6m，y方向的最大位移解析解(uy)max=2.446 4×10-5m，法向应力最大值解析解(σxx)max=5×104Pa.

3.3 有限元分析计算

采用ansys workbench对悬臂梁模型进行有限元分析，网格划分结果如图5 所示，计算结果如图6 所示.

图5 网格划分结果Fig.5 Meshing result

(a) x方向位移变化云图(b) y方向位移变化云图(c) 法向应力分布云图图6 有限元分析结果Fig.6 Analysis results of finite element

3.4 等几何分析计算

计算u方向的节点向量，其中，n=5,p=3,m=9.

利用式(5)和式(6)计算d：

|Q1-Q0|=6,|Q2-Q1|=6,

|Q3-Q2|=6,|Q4-Q3|=6,

|Q5-Q4|=6,

d=30.

利用式(8)计算uj+p：

u0=…=u3=0,u6=…=u9=1，

计算v方向的节点向量，其中，n=2,p=2,m=5，直接根据式(8)可知

v0=…=v2=0,v3=…=v5=1，

所以，v=[0 0 0 1 1 1.

采用matlab编写计算过程. 首先建立NURBS基础数据库，如表1 所示.

通过调用函数生成控制点网格和节点向量网格(其中经过3次插值细分)，如图7 和图8 所示.

通过matlab编程初始化边界，计算单元刚度矩阵并组装，再通过Gauss积分计算节点力，边界条件处理调用求解，最后的后处理结果如图9 所示.

表1 NURBS基函数Tab.1 NURBS basis function

图7 控制点网格

图8 节点网格Fig.8 Node grid

(a) x方向位移变化云图

(b) y方向位移变化云图

(c) 法向应力分布云图图9 等几何分析结果Fig.9 Isogeometric analysis results

3.5 计算结果对比

对比有限元分析(FEA)和等几何分析(IGA)计算的x方向和y方向的位移变化以及法向应力(见图6 和图9)可知，两种方法所得位移变化和法向应力分布云图基本一致，且最大值点位置相同，最大值与解析解误差如表2 所示.

表2 FEA与IGA法的计算结果对比Tab.2 Comparison of calculation results in FEA and IGA

在相同的计算条件下，由matlab编程所产生的节点数量和网格数量等于或者大于由有限元分析生成的节点数量和网格数量. 为了保证两种方法运算时间对比的公平性，采用控制变量，选择与解析解精度相同的有限元解，保证了求解结果和解析解的一致性，从而使获得的运算时间更具有代表性，同时，为了减小计算机硬件资源对运算效率分析的影响，特别是AVX指令集在CPU频率不同时的运算效率差距很大，将频率固定在3.0 GHz，内存频率固定在3 200 MHz，以减少自动睿频对分析结果的影响. 分析时间定义为：有限元分析时间从网格划分结束后开始到后处理结束(人为设置边界条件和后处理方案不计入分析时间)；
等几何分析时间指全部分析、计算、后处理时间. 总时长定义为：网格划分时间和分析时间的总和. 每个算例中，计算所耗费的时间都取10次的平均值. 运算时间统计结果如表3 所示.

表3 FEA与IGA法的运算时间对比Tab.3 Comparison of calculation time in FEA and IGA

由表2 可知，本文改进后的等几何分析计算的x方向位移变化、y方向位移变化以及法向应力变化比有限元法计算的精度更高. 由表3 可知，在保证计算精度的前提下，本文等几何分析计算的总时长比有限元分析计算的总时长缩短了47.72%.

本文利用平均值法对等几何分析的节点向量求解进行改进，建立了改进等几何分析处理平面静力学问题的结构化流程. 该流程用平均值法代替原等距分布法，建立NURBS基础信息库，然后通过等几何建模、单元细分、施加边界条件、后处理计算四个步骤进行求解，并将其应用于分析经典静力学问题. 通过等几何分析法求解悬臂梁问题得到的x方向、y方向位移变化云图与有限元分析一致，且最大值与解析解相比的误差分别为0.84%和0.06%，法向应力最大值的误差比有限元的误差减小了约5%，验证了本文方法的可行性；
在保持计算精度的情况下，改进后的等几何分析的运算总时长比有限元分析缩短了47.72%，证明了本文方法的高效性.

猜你喜欢样条控制点平面对流-扩散方程数值解的四次B样条方法图学学报(2020年5期)2020-11-13立体几何基础训练A卷参考答案中学生数理化·高三版(2019年1期)2019-07-03三次参数样条在机床高速高精加工中的应用制造技术与机床(2017年7期)2018-01-19NFFD控制点分布对气动外形优化的影响北京航空航天大学学报(2017年4期)2017-11-23三次样条和二次删除相辅助的WASD神经网络与日本人口预测软件(2017年6期)2017-09-23基于样条函数的高精度电子秤设计计算机测量与控制(2017年6期)2017-07-01基于风险管理下的项目建设内部控制点思考中国工程咨询(2017年12期)2017-01-31参考答案试题与研究·高考数学(2016年1期)2016-10-13关于有限域上的平面映射肇庆学院学报(2016年5期)2016-03-11SDCORS在基础地理信息控制点补测中的应用全球定位系统(2015年4期)2015-02-28

推荐访问:静力学 平均值 几何