局部近似特别解的无网格方法

局部近似特别解的无网格方法

在传统的基于网格的微分方程数值方法（如有限差分法（fdm）和有限元法（fem）中，传统的基于网格的数值法广泛应用于计算水流和固体力学领域。但这些方法仍然存在诸多不足之处,因为其在很大程度上取决于网格划分的质量,当划分的网格出现畸变时会严重影响数值解的精度甚至导致计算失败,从而使得这些方法在许多问题上的应用受到了限制。于是无网格方法应运而生,并且在求解偏微分方程方面得到很大的关注。无网格方法包括:Lucy和Gingold提出的光滑粒子流体动力学方法(SmoothedParticleHydrodynamics,SPH),基本解方法(MFS),局部近似特别解方法(LMAPS)等。在众多的无网格方法中基于径向基函数(RadialBasisFunctions,RBFs)的无网格方法备受瞩目。RBFs是一类以场点和源点之间的欧式距离为自变量的函数,它具有形式简单,各项同性等优点,已被成功应用于多变量插值。但基于全局性质的RBFs的无网格方法产生的系数矩阵常常是满秩,有时甚至是病态的,适于解决有限插值点的问题,但应用到解决大规模散乱点的数值计算问题时则受到限制。为了克服上述问题,文中采用局部性质的近似特别解方法解决大规模问题,研究表明效果显著。该方法只需创建局部区域并利用该区域的点构造局部低阶矩阵,然后将局部形式推广至全局形式形成大规模稀疏矩阵,最后通过求解这个稀疏线性方程组就可以得到偏微分方程的近似解,从而可以成功避免上述问题。1局部化的map方法首先考虑如下的椭圆形偏微分方程:其中,Δ为二阶线性Laplace算子;a(x),b(x),c(x),f(x)和g(x)为给定的函数;B为边界算子。MAPS的主要思想是将式(1)转化为Poisson型方程,即其中,H(x,u,ux,uy)=-a(x)ux-b(x)uy-c(x)·u+f(x),并且H可以用近似表示运用中常取MQ函数作为基函数,其表达形式为其中,c为形参,决定基函数的形状。当c很小时,MQ函数成锥形形状;当c逐渐变大时,函数图像则越来越平滑。由此可以看出,形参c对于近似式的精确性至关重要。影响c的取值因素主要有3个,即支撑域内节点的数量、支撑域内的节点分布以及支撑域的大小。一般情况下,c越大,近似误差就越小。应用MQ函数插值近似会出现这样的情况:为提高近似精度增大c值时得到的系数矩阵往往会变得越来越病态,因而舍入误差会变大,增加了解的不稳定性。此外,如何适当选取c值并不是一件容易的事情,这也成为学界研究讨论的热点问题。根据式(5)和式(6)推得由式(3)、式(4)和式(7)可知式中:重组式(8),有其中并且式(2)可写为下面具体介绍局部化的MAPS方法。由于Ωi={xki}kni=1中的每一个点对式(11)都成立,于是有在式(12)中右端的方阵可以用Φni表示,并且Φni是非奇异的,因此是可逆阵,所以得到的式(12)中未知系数为其中因此,在式(11)中的^u(xi)可以表示为ni个节点函数值的线性组合。其中类似于式(15)有对于边界条件,在局部区域上,有于是在全局区域上有所以由式(17)和式(19)可以构造一个N维稀疏线性方程组,通过求解该方程组可以得到求解该区域上所有结点上的近似解。2局部近似方法利用LMAPS进行数值计算时,找到区域中每个计算结点相邻近的ni个结点是很重要的,当插值点为很庞大的数值时,有效的搜索算法非常重要。在这些算法中,kd-tree算法是一种k维数据空间的有效搜索算法。为了验证文中数值精度,分别定义均方根误差(RMSE),最大绝对误差(MAE),最大相对误差(RAE)如下:式中:N为测试点的数目;^uk为近似解;uk为精确解。算例1考虑Dirichlet边界型条件的Poisson方程:其中,Ω∪ue785Ω为如图1所示的一个不规则区域。由图1可以看到,内部和边界结点都是均匀分布的,在计算时取14408个内部结点和300个边界结点。该二维Poisson方程的精确解为u(x,y)=e-2x+3y。图2为分别利用LMAPS方法和传统的LDQ方法得到的绝对误差。由图2(a)和图2(b)比较可以看出,LMAPS方法得到的误差比LDQ方法更小,因此,LMAPS方法略优于LDQ方法。同时还可以看出,误差在不规则区域内的大部分情况都是下降平稳的,除了在边界周围的一部分点有一定的跳跃。但作为局部方法,它们与全局的径向基函数方法比较,计算效率还是比较显著的。算例2采用局部近似特别解方法求解方腔驱动流问题,以检验这种方法的精度和准确性。不可压流体是由如下原始变量形式的Navier-Stokes方程控制。动量方程:连续方程:图3为Re=50和Re=500时的方腔驱动流的流线图。当Re很小时,靠近腔体上壁有一个呈顺时针方向的漩涡;随着Re的增大,方腔内部受到流体流动的影响也会变大,从而使得靠近壁面的速度梯度越来越大且腔内主漩涡逐渐向腔体右下方移动。当Re=50腔体右下角和左下角先后次生出一个小涡;随着Re继续增大,这两个次生的小涡会越来越大,而主漩涡则逐渐向腔内中心位置移动。图4为不同雷诺数下方腔垂直中心线上速度u和水平中心线上速度v的分布情况。根据文献将得到的数值解与Ghia和Sanyasiraju的结果进行比较,得到的解和参考解具有相似的精度。如果增加区域内结点分布的密度,局部近似特别解方法的精度还会有一定的提高。由图3中可以看出,涡流开始在方腔中心轴线的上半部分且处于一个稳定的状态,随着雷诺数的增大涡流逐渐向方腔的中心移动。由图4可以看出,雷诺数从50到500,在薄边界层的附近,u,v的速度梯度在不断地增加。3大力采用