二次有限元格式高阶数值模拟矢量型弹性方程stress边值问题

刘 雪, 丁 晓, 江 山

(南通大学理学院, 江苏 南通 226019)

弹性力学主要研究弹性体受外界作用、环境变化或边界约束等因素影响所导致的应力、应变和位移变化情况,被广泛应用于道路桥梁、水利工程和航空航天等领域. 近年来, 处理弹性力学方程的数值方法主要有有限元法[1-3]、混合有限元法[4-6]、弱有限元法[7-9]和多尺度方法[10-12]等, 其中有限元法因具有诸多优点而备受关注: 可针对不同方程建立不同的变分形式; 等参单元与实际单元之间易于实现仿射变换;能通过独立构造分片基函数有效耦合局部与全局的数据信息; 数学理论完备等. 然而, 已有的工作主要采用线性基函数构造多项式, 且多为处理弹性力学方程的本征边界类型,抑或是格式较为复杂,计算成本偏高. 本文拟采用有限元法的高次基函数格式, 精确模拟带有复杂stress边界条件的矢量型弹性方程,通过法向量和切向量有效刻画对应的分量, 在更精细的数据结构下组装成总刚度矩阵并求解代数方程组, 以期利用有限的计算资源得到高阶精度且稳定收敛的数值模拟.

1 弹性方程及其变分

考虑二维弹性力学方程

(1)

(2)

(3)

2 有限元基函数的构建

B(uh,vh)=(f,vh), ∀vh∈Uh0×Uh0,

(4)

其中左端的双线性形式

(5)

右端项(f,vh)需要根据具体的边界条件进行相应的调整, 这里Uh0表示边界为零的紧支集函数空间.

分别选用Lagrange型线性基函数和二次基函数构造有限维泛函空间.将区域离散化剖分形成三角形网格，然后以每个三角形单元的3个顶点作为节点构造局部线性基函数.在此基础上, 再增设3条边的中点作为单元新节点,进而构造局部二次基函数.

1) 线性基函数的构造.选取等参三角形单元的3个顶点Α1=(0,0),A2=(1,0),A3=(0,1),A1,A2,A3须满足局部基函数ψj的构造条件:

2) 二次基函数的构造.除选取等参三角形单元的3个顶点Α1=(0,0),A2=(1,0),A3=(0,1)之外, 再增加3条边的中点A4=(1/2,0),A5=(1/2,1/2),A6=(0,1/2).各节点均须满足局部基函数ψj的构造条件:

3 有限元法处理stress边界条件

(6)

进一步地,变分形式(3)为

(7)

(8)

根据式(7)可知, 在处理stress边值时须先对其右端标量积项进行细节处理.对于局部边界Γs⊂∂Ω, 设定二维问题边界Γs的4条边的法向量依次为nt=(0,1)T,nb=(0,-1)T,nl=(-1,0)T,nr=(1,0)T, 切向量依次为τt=(-1,0)T,τb=(1,0)T,τl=(0,-1)T,τr=(0,1)T. 将其代入具体边界条件, 可得出pn和pτ.由式(8)可知, 右端的载荷向量还存在切线方向和法线方向的stress边值.故在法线方向和切线方向进行曲线积分, 分别对应添加到f的两个分量f1,f2上, 从而组装为新的总载荷向量, 即全局区域的stress边界总载荷向量.

4 误差分析

引理1[13]若问题(1)的真解u和右端f具备必要的光滑性, 则线性有限元解uh的误差估计为

‖u-uh‖∞≤Ch2‖f‖0,

(9)

‖u-uh‖0≤Ch2‖f‖0,

(10)

其中C为与单元步长h无关的常数.

引理2在引理1条件下, 线性有限元解uh的H1半范数误差估计为

|u-uh|1≤Ch‖f‖0.

(11)

定理1在引理1条件下, 二次有限元解uh的误差估计如下:

‖u-uh‖∞≤Ch3‖f‖0,

(12)

‖u-uh‖0≤Ch3‖f‖0,

(13)

|u-uh|1≤Ch2‖f‖0.

(14)

证明 限于篇幅, 这里仅针对L2误差范数证明式(13).二次有限元格式取k=2的分片二次多项式, 解向量u的分量ui(i=1,2)均有‖ui-uih‖0≤Ch|ui-uih|1, 而|ui-uih|1≤Chk|ui|1, 则‖ui-uih‖0≤Chk+1|ui|1.进一步根据引理2, 式(13)得证.

5 数值实验

针对带stress边界的矢量型弹性方程, 通过MATLAB编程运行算例的实际结果, 利用误差图示与数值分析,分别采用线性有限元格式和二次有限元格式验证本文有限元方法的模拟效果.

当二维问题各方向取剖分数为26时, 其x方向等距步长h=2-6, 分别应用线性有限元格式和二次有限元格式计算得到解向量的两个分量u1,u2的绝对误差如图1～2所示.由图1～2可见: 线性有限元格式下u1,u2的最大误差上限为4×10-5, 而二次有限元格式下的最大误差上限改进为4×10-7,表明本文方法的精度更高;相较线性有限元格式,二次有限元格式的误差的最大值、宽度和高度都明显减小,故二次有限元格式求解带stress边界的矢量型弹性方程时精确性更高且收敛效果更优越.

图1 h=2-6时线性有限元格式下分量u1(a), u2(b)的绝对误差Fig.1 Absolute errors of components u1(a), u2(b) in linear finite element scheme, with mesh size h=2-6

图2 h=2-6时二次有限元格式下分量u1(a), u2(b)的绝对误差Fig.2 Absolute errors of u1(a), u2(b) in quadratic finite element scheme, with mesh size h=2-6

为了突出对比线性有限元格式与二次有限元格式的具体精度和收敛效果, 对初始粗网格步长h=2-2进行折半递减的剖分加密.在各自粗细网格分别计算矢量型真解与有限元解向量之间误差的L∞范数、L2范数和H1半范数, 进而计算对应的收敛阶,结果如表1～2所示.由表1～2可见: 在相同的步长剖分下,二次有限元格式较线性有限元格式的精度高达103倍.随着区域离散三角形网格的加密,线性有限元解与真解之间误差的L∞范数和L2范数为二阶收敛,H1半范数为一阶收敛;而二次有限元解与真解之间误差的L∞范数和L2范数均实现了三阶收敛,H1半范数为二阶收敛,与定理1的误差分析结果完全一致.这充分验证了二次有限元格式的误差范数值显著降低,收敛阶和收敛速度显著提升.

6 结语

本文基于有限维空间逼近无限维空间的数值思想,利用有限元法处理具有stress边值的二维矢量型弹性力学方程. 有限元格式在每个三角形单元分别构造线性基函数和二次基函数, 对应形成稀疏程度和数据结构不同的代数方程组, 利用法向量和切向量刻画stress边值对应的函数细节. 理论证明与数值实验验证了本文方法的高阶稳定收敛,充分展现了高次有限元计算格式的高阶模拟优势,为有效求解弹性边值问题提供了高效数值逼近的参考方案.

表1 各网格步长下线性有限元格式的误差范数值及其收敛阶

表2 各网格步长下二次有限元格式的误差范数值及其收敛阶