刘 雪,丁 晓,王晓莹,江 山
(南通大学 理学院, 江苏 南通 226019)
弹性力学在科学工程中的应用极为广泛,比如道路桥梁、房屋建筑、水利建设、航空航天等领域,其有效求解已成为解决实际工程问题必不可少的手段[1-4]。目前,主要研究方法有物理实验法、有限元法、有限差分法、边界元法、无网格法等。近年来,Tarasov[5]研究了梯度弹性力学方程,结合差分法和物理法进行求解;王芬玲等[6]利用双线性元与零阶R-T 元,面向拟线性粘弹性方程建立了混合元数值格式;孙美玲[7]使用有限元法结合自适应分层网格有效逼近了奇异摄动的边界层问题,得到了一致收敛的高精度结果;Chung等[8]研究广义多尺度有限元法处理静态线性玻尔兹曼方程;文献[9]应用波函数法求解了弹性动力学的二维各向异性的时间调和方程;文献[10]针对线性粘弹性对称模型,得到了混合有限元法的具体实现。
本文将应用有限元法处理弹性力学方程,利用有限维逼近无限维的数学思想,分别采用线性元基函数与二次元基函数进行数值逼近。随着二维离散网格的剖分加密,不断得到数值精度更好、收敛阶更高的实践结果。
考虑二维方程
下面将原方程转化为变分形式。使用虚功原理,方程(1)两端同时乘以向量型检验函数Green 公式和分部积分,依据本征边界条件限制一阶可导且平方可积的紧支集函数空间。可知
故原方程(1)对应的变分形式为:
寻找u→∈H1(Ω)×H1(Ω),使得
数值计算方法旨在通过有限维空间近似无限维空间,保证数值解有效近似真解,进而考虑其数值稳定等性质。
有限元法利用分片多项式构造基函数,进而由基函数张成相应的有限维泛函空间,比如选择一组有限维基函数 φ1,φ2,…φn,记有限元空间程(1)是张量方程,其目标解向量u→ =(u1,u2)T具有二维分量,用数值解向量u→h=(u1h,u2h)T表示近似,应的变分形式为寻找u→h∈Uh×Uh,使得
其中,Uh0表示边界上为零的紧支集函数空间;即为张量形式的有限元近似解,需通过有效数值格式求解其各自的分量u1h、u2h,用以逼近真解u→的各自分量 u1、u2。
对求解区域Ω 进行网格剖分,选用三角形单元做区域离散化处理。通过仿射变换来建立任意三角形单元和等参单元之间的对应关系,从而构造局部基函数ψ。此处先构造线性元局部基函数,选取三角形的三个顶点为对应坐标,此时对应的局部基函数分别为ψ1= 1 - x - y,ψ2= x,ψ3=y,于是每个单元有三个自由度。这样在某个
最后,通过全局节点与局部节点的对应关系,形成总刚度矩阵A 和总载荷向求解代数
构造数据结构更复杂、逼近效果更优越的二次元基函数,进而由局部单元映射到全局单元,与线性元基函数相比,会形成稠密程度稍高的总刚度矩阵A,但其依然为稀疏矩阵。
构造二次元局部基函数时,在等参三角形单元上除了同样选取三角形的三个顶点外,再选三条边的中点作为对应坐标,从而形成每个单元有六个自由度的二次元基函数。通过运算,局部基函数分别为:
通过编写Matlab 程序和算例实践来验证文中方法的正确性,给出数值精度和收敛阶数,更直观详实地比较线性元格式与二次元格式的优劣。定义真解u→和有限元解u→h之间误差的 L∞范数、L2范数与H1半范数分别为
图1 的左右两子图直观画出了线性有限元解逼近真解的绝对误差三维图,解向量的两个分量最大误差都是2×10-5数量级,可见数值逼近真实可靠。为了进行公平比较,图2 进一步显示了在同样剖分数N = 64 条件下,二次有限元解有更加精确的逼近效果,解向量的分量上限均继续降为4×10-9数量级,表明二次元较之线性元的改进效果直观且显著。
图1 用线性有限元格式求真解u1、u2 对应的绝对误差(取剖分数N=64)
图2 用二次有限元格式求真解u1、u2 对应的绝对误差(取剖分数N=64)
为更好地展现有限元法采用线性基函数与二次基函数的联系和区别,在单方向剖分数以2n逐渐加密离散化网格的基础上,观察该二维问题。依式(5)-(7)计算随着一致网格剖分数的不断倍增,真解与有限元解之间误差的L∞范数、L2范数以及H1半范数的具体数值,给出随之产生的上下粗细两层网格的收敛趋势,结果如表1、表2。表1显示,在选取线性元基函数时,L∞范数、L2范数都有二阶收敛,H1半范数有一阶收敛;表2 显示,选取二次元基函数时,其误差范数值较之表1 显著降低,L∞范数、L2范数都实现了更高的三阶收敛,H1半范数实现了更高的二阶收敛。横向与纵向对比都表明,二次有限元计算格式具有更好的计算精度、更高的收敛阶数。
表1 线性有限元格式的误差范数
表2 二次有限元格式的误差范数
综上所述,通过本文的研究与实践,充分展现了有限元法求解弹性力学方程解向量的可靠效果,而分别采用线性基函数、二次元函数,可进一步展现了高次元计算格式得到的高精度、高收敛的一致稳定数值结果。