多重网格法在非结构网格中的应用

多重网格法在非结构网格中的应用（精选8篇）

篇1：多重网格法在非结构网格中的应用

多重网格法在非结构网格中的应用

将多重网格法运用于非结构网格.网格是通过聚合法得到的,网格之间是相互关联的.方程的`求解采用Jamson的有限体积法.给出了二维、三维情况的数值算例.

作 者：刘学强 伍贻兆 夏健  作者单位：南京航空航天大学空气动力学系.江苏,南京,210016 刊 名：计算物理  ISTIC EI PKU英文刊名：CHINESE JOURNAL OF COMPUTATIONAL PHYSICS 年，卷(期)： 19(4) 分类号：V211.3 关键词：多重网格   有限体积法   非结构网格  

篇2：多重网格法在非结构网格中的应用

随着计算流体力学的发展,多重网格在求解欧拉方程、N-S方程过程中得到了广泛应用.将这种方法应用于高超声速进气道的数值模拟中,通过具体的算例,验证了该方法的.有效性.

作 者：尹凤娥 YIN Feng-e  作者单位：中南林业科技大学,土木工程与力学学院,湖南,长沙,410004 刊 名：航空计算技术  ISTIC英文刊名：AERONAUTICAL COMPUTING TECHNIQUE 年，卷(期)： 40(2) 分类号：V211 关键词：多重网格   高超声速   有限体积法  

篇3：多重网格法在非结构网格中的应用

无网格法是近二十年来才发展起来的数值方法[1],该方法可以部分或彻底消除网格,而采用基于点的近似,则可以完全抛弃网格重构。目前应用最广泛的无网格法是无单元伽辽金方法[2],其特点在于:1)采用移动最小二乘法近似位移函数。2)从能量泛函的弱变分形式得到控制方程。3)采用拉格朗日乘子法[3]满足本征边界条件。

2 无网格法基本原理

从数学理论上讲,最小二乘法描述的是函数逼近问题,即最小二乘逼近。假设待求函数u(x)在求解域Ψ中的N个节点xI(I=1,2,…,N)处的函数值uI=u(xI)是已知的,在域Ψ内构造待求函数u(x)的全局近似函数uh(x)。待求函数u(x)在计算点x(在伽辽金无网格法中为高斯点)的邻域Ψx内可以局部近似为:

其中,为计算点邻域Ψx内各点的空间坐标,;为基函数;m为基函数的个数;a(x)=[a1(x),a2(x),…,am(x)]T;ai(x)为待定系数。

在移动最小二乘近似(MLS)中,系数ai(x)的选取使得近似函数uh(x,)在计算点x的邻域Ψx内是待求函数u(x)在某种最小二乘意义下的最佳近似。

设计算点x的邻域Ψx包括N个节点,近似函数uh(x,)在这些节点处的误差的加权平方和为:

令J取最小值,即:

其中形函数N()为:

近似函数uh(x,)是待求函数u(x)在计算点x邻域Ψx内的加权最小二乘意义下的局部最佳近似。

3 无网格法的弹塑性理论

3.1 屈服条件

Drucker-Prager屈服条件的表达式为:

其中,α,k均为正常数。

平面应变问题:

其中,c为粘聚力;φ为内摩擦角。

3.2 加载和卸载准则

理想弹塑性材料[4]不发生强化,其屈服面的大小和形状不随内变量的发展而变化,因而只有F(σij)=0屈服面,当然也可认为后继屈服面与初始屈服面重合。

在这样的条件下,F(σij)<0表示材料处于弹性状态;F(σij)=0,同时发生应力增量dσ时,有两种不同反应:

1)应力点保持在屈服面上,这时有新的塑性变形增量dεP产生,这种情况称为加载;2)若应力点由屈服面退回屈服面内,则称为卸载,其表达式如下:

3.3 弹塑性本构矩阵的一般表达式

增量塑性理论要对以下三方面作出基本假定:1)屈服条件;2)流动法则;3)硬化法则。设屈服条件用下式表示:

其中,σij为应力状态;K为硬化函数。

在增量理论中,材料达屈服以后把应变增量分为弹性增量和塑性增量两部分,即:

其中,弹性应变部分与应力增量之间的关系仍服从胡克定律,即:

其中,[D]为弹性矩阵。

应用流动法则得:

此即为增量理论的弹塑性矩阵通式。

4 无网格法求解土质边坡稳定

算例验证:均质土坡,坡高H=20m,土容重r=25kN/m 3,粘聚力c=42kPa,内摩擦角φ=17°,坡角β=30°,35°,40°,45°,50°。

如图1所示,按照平面应变建立模型(β=45°),边界条件为左右两侧水平约束,下部固定,上部为自由边界,采用关联流动法则。

图2 为该边坡有限元(ANSYS)计算收敛(泊松比υ=0.3)时的塑性区分布,显然此时边坡的塑性区已经贯通,但是边坡并没有产生很大的塑性应变和位移,而是处于稳定状态。图3为无网格法计算的滑裂面形状。从图2,图3中可以看出,有限元法与无网格法所计算出的滑裂面的位置是非常接近的。

5 结语

无网格法能够有效地对土质边坡稳定问题进行有效的弹塑性分析。与有限元法相比,无网格法不需要划分单元,只需布置节点,因此特别适合于各种裂纹问题以及弹塑性问题的分析,如本文讨论的土质边坡稳定问题。本文计算了土质边坡稳定的安全系数,数值结果表明了无网格法求解弹塑性问题的可行性。

参考文献

[1]张雄,刘岩.无网格法[M].北京:清华大学出版社,2004.

[2]周维垣,宼晓东.无单元法及其在岩土工程中的应用[J].岩土工程学报,1998,20(1):5-9.

[3]刘天祥,刘更.二维接触问题的无网格伽辽金—有限元耦合方法[J].西北工业大学学报,2003,21(4):499-502.

篇4：多重网格法在非结构网格中的应用

关键词：多重网格；温纳；限制算子；插值算子

1.多重网格法的发展现状

多重网格法简称MG，是数学家Fedorrenko于20世纪60年代最早提出的，直到Brandt关于“边值问题多重网格适应解”的文章发表后，多重网格法才开始得到广泛的关注和全面研究[1]。多重网格法是目前近乎最优的算法之一，有其固有的特点：收敛速度跟网格的大小无关，而且这一特性同样适应于求解复杂问题。多重网格法经过近30年的研究和发展，有了完整的理论体系，算法趋于成熟，在流体力学，计算电磁学，波动方程等领域都有所应用，是当今数值计算领域常用的算法之一。近二十年也逐渐的应用于地球物理领域中。

2.多重网格法的基本思想及理论公式

多重网格法其主要思想是在不同的网格层中，可以使不同频率的误差分量因网格间距的不同逐渐被衰减掉[2]。首先把以当前网格为基础的低频分量限制到比较粗的网格上，这样就会使原本看起来比较光滑的的低频分量在粗网格上又重新跳动振荡起来，然后根据松弛迭代的方法对其进行处理，这样会使在粗网格上误差又变得光滑起来，继续将其限制到更粗的网格上进行光滑迭代，重复以上过程直至最粗限度的网格。到最粗的网格结束后将其每层网格迭代光滑后的误差经线性插值后再逐级返回至细一层网格，将这些返回后的误差对方程解进行校正，就可以得到更为精确的解，至此一次多重网格的循环结束[3]。

以2层网格为例进行理论推导：

假设存在一个均匀剖分的两重网格如图，细网格的步长是h，粗网格的步长是H，H=2h，用线性方程表示剖分最细的网格层：

图1 二维模型均匀剖分图

Ahuh=fh（1）

Ah是系数矩阵，fh是右端项，uh是方程的精确解。

在细网格层上进行1次松弛迭代，uh1是迭代过后的近似解，uh1和uh间高频分量因迭代而衰减。

uh-uh1=eh1是误差向量，rh1就是残差向量

Ah（uh-uh1）=Aheh1=fh-Ahuh1=rh1（2）

将（2）式近似到粗网格上，AHeH=rH1（3）

由细到粗需要加入限制算子RHh（线性算子），使rH1=RHhrh1（4）

求解（3）得到eH1，加入插值算子PhH（线性算子），

得到eh1的近似解e1h，则e1h=PhHeH1（5）

可以方便的计算出公式（4）和（5）。最后就可以得出修正后的近似解h1

h1=uh1+e1h（6）

这是2层V循环的理论公式，亦可以看成任意相邻网格层间的转换公式。

在上述求解过程中，限制和插值算子是其重要组成，常用的是线性插值及完全加权限制，因是均匀剖分，故其常用的9点限制算子和插值算子分别是：

RHh=1161811618141811618116，PhH=121242121（7）

3.MG法三维响应模拟

将多重网格法应用于以有限差分为基础的高密度正演中，以温纳装置对典型的地下异常体进行测量，电流都是I=1A，电极距为5米，测线总长度100米，则测区范围一共有21个电极。网格剖分为49×49×31，剖分模式为不等距剖分，中间细网格，边界剖分为指数增长的网格。构建三维模型

3.1不同测线方向的异常响应

图2 不同测线分布图

为了突显正演的响应，在用温纳测量时，从不同的倾角进行测量，如图2θ是测线与X轴的夹角，（θ0=0；θ1=30；θ2=45；θ3=60；θ4=90），在用温纳装置测量时，会得到63个测点数据，测得地下深度42.5米；在用偶极-偶极测量时，会得到171个测点数据，测得地下深度47.5米。以低阻异常为例，异常体的电阻率为1Ω·m，围岩电阻率为100Ω·m，异常体埋深为Z方向3到10米，X方向为-10米到10米。

当地下异常体是一个沿Y方向延伸的二度体，通过温纳装置从不同角度的到的视电阻率图（图3）。

图3 Y=100：温纳不同角度视电阻率图

利用温纳装置测量时，当θ0=0°时测线垂直横穿过低阻异常体，此时可以得到一个范围比较小可近似看成倒梯形的低阻异常范围，实际范围X方向比较准确在模型范围内，Z方向向下延伸比较多；当θ0=30°、45°、60°时同样可以得到异常区域范围X方向逐渐变大，直至90°时通过图e可以看到，视电阻率图所呈现的是一个地下均匀平行的两层介质。

3.2异常体埋深的响应

前面得到了测线穿过异常体的各种情况，本节将研究当模型参数不变时，测线不穿过异常体，并与异常体垂直距离不等时的响应特征。温纳的装置参数都不变，当异常体垂直距离测线Y=0、15、35米时如图4。

经过两种装置测量后，得到图5，如图所示，随着与测线垂直距离的不断加大，正演响应均有一定变化。当Y=0m时，相当于垂直横穿过异常体，与5.1中θ=0时情况相同；当Y=15m时，异常区域均有向下的趋势，但是X方向范围变化不大；当Y=90m时，异常区域都为位于图片的下方，几乎脱离视电阻率图。总结来说，当异常体与侧向越近时，响应特征越好，随着距离的增加，异常响应会逐渐向下，最后移出测量范围。

图4 不同埋深时示意图

图5 不同埋深时温纳装置视电阻率图

结论

在应用多重网格法进行迭代计算的基础上，用温纳装置对地下典型的异常体模型进行正演模拟，通过对正演响应的分析发现，地下存在二度体和异常体Y的长度变短时的情况下，从不同角度进行测量时响应的结果差别很大，只有当测线垂直穿过异常体时，响应结果最好；当异常体与测线距离不同时，距离越近，正演响应越接近真实模型。以上两种正演响应均能较好的反应出不同条件下的模型特征，相信多重网格法未来在物探领域应该有更好的发展。

参考文献：

[1]Alcouffe R.E.，Brandt A.，Dendy J.E.，and Painter J.W..The multi-grid method for the diffusion equation with strongly

[2]乔中林.直流电法三维正演多重网格算法研究[D].北京：中国地质大学，2007.

[3]鲁晶津.地球电磁三维数值模拟的多重网格方法及其应用研究[D].合肥：中国科学技术大学，2010.

篇5：多重网格法在非结构网格中的应用

以三维非结构网格的显式有限体积法为基础,采用了一种TVD方法求解三维Euler方程,使用Roe通量来计算网格单元边界处的守恒量通量.为了验证方法的.可行性,用该方法模拟三维爆炸问题,得出的结论是我们的方法可行,稳定且有效.

作 者：高二 宋松和 GAO Er SONG Song-he  作者单位：国防科技大学,湖南,长沙,410073 刊 名：航空计算技术  ISTIC英文刊名：AERONAUTICAL COMPUTING TECHNIQUE 年，卷(期)： 38(5) 分类号：V211 关键词：TVD方法   双曲型守恒律方程组   三维非结构网格   有限体积法  

篇6：多重网格法在非结构网格中的应用

由于变后掠翼可以更好地解决高、低速性能要求的矛盾, 是现代飞机设计的一个重要途径, 所以选此作为研究对象, 进行初步探索性研究。

机翼变后掠流场的求解包含运动边界的非定常复杂流场, 是计算流体力学中的难点。非结构动网格技术 (DUT) 是在非结构网格框架下求解此类问题的关键。目前国内外发展的可靠性较高的网格变形方法主要有弹簧近似方法[1]和弹性体方法[2], 弹性体方法的网格变形能力强且网格质量较高, 但计算工作量大;弹簧近似方法则比较简单、高效, 但网格变形能力不如弹性体方法。

通过对弹簧近似方法的改进, 使弹簧近似方法在保持原有计算效率的前提下, 网格变形能力和网格质量得到大大的提高。最后, 将发展后的非结构动网格生成方法与ALE形式的流场解算器耦合起来, 求解了包含运动边界的、模拟上述机翼运动状态下的亚, 跨音速绕流, 并对结果进行了分析。

1 弹簧近似方法

根据Blom[3]的定义, 弹簧近似方法有两种描述方法, 一种是弹簧的平衡长度为零, 称为顶点弹簧 (Vertex springs) ;另一种是弹簧的平衡长度等于边界未运动和变形前边的原长, 称为棱边弹簧 (Segment springs) 。

传统的顶点弹簧近似方法一般用于网格优化, 弹簧倔强系数取为常数, 网格点移动后的新位置坐标相当于周围点坐标的平均。现对这一经典方法进行了修改, 使其更适于实现网格变形运动。下面给出修改后的顶点弹簧近似方法。对于顶点弹簧, 节点i, j间的弹簧张力可以写为:

$\overrightarrow{F}{}_{ij}={\rm K}{}_{ij}(\overrightarrow{r}{}_j-\overrightarrow{r}{}_i)(1)$

(1) 式中Kij为连接节点i, j的弹簧的倔强系数;$\overrightarrow{r}$i, $\overrightarrow{r}$j分别为节点i, j的位置矢量。对任意的一个节点i, 它的合力方程为

$\overrightarrow{F}{}_i={\displaystyle \sum _{j=1}^{{\rm N}{}_i}{\rm K}}{}_{ij}(\overrightarrow{r}{}_j-\overrightarrow{r}{}_i)(2)$

(2) 式中Ni为计算区域内所有与节点i相连的非结构网格的节点总数。不同于网格优化的情况, 修改后的顶点弹簧方法, 令计算区域内所用网格节点的受力都为边界没有发生运动和变形时所受的合力。对计算区域内所有的m个非结构网格的节点, 列上述合力方程, 得到初始状态下的该线性系统的矩阵表达式

$\begin{array}{l}\left[\begin{array}{cccccc}-{\displaystyle \sum _{j=1}^{{\rm N}{}_i}{\rm K}}{}_{1j}& \cdots & {\rm K}{}_{1j}& \cdots & {\rm K}{}_{1{\rm N}}& \\ ⋮& \ddots & \cdots & \ddots & ⋮& \\ {\rm K}{}_{i1}& \cdots & -{\displaystyle \sum _{j=1}^{{\rm N}{}_i}{\rm K}}{}_{ij}& \cdots & {\rm K}{}_{i{\rm N}}& \\ ⋮& \ddots & \cdots & \ddots & ⋮& \\ {\rm K}{}_{m1}& \cdots & {\rm K}{}_{mj}& \cdots & -{\displaystyle \sum _{j=1}^{{\rm N}{}_i}{\rm K}}{}_{mj}& \end{array}\right]\times \\ \left[\begin{array}{c}\overrightarrow{r}{}_1\\ ⋮\\ \overrightarrow{r}{}_i\\ ⋮\\ \overrightarrow{r}{}_{{\rm N}}\end{array}\right]=\left[\begin{array}{c}\overrightarrow{F}{}_1\\ ⋮\\ \overrightarrow{F}{}_i\\ ⋮\\ \overrightarrow{F}{}_m\end{array}\right](3)\end{array}$

(3) 式中Ni为包括边界点的网格节点数。当边界运动或变形后, 通过改变边界的位置矢量和适当的整理, (3) 式就得到一个m阶的线性方程组。由于该线性方程组的刚度矩阵对称正定, 则求解的Gauss-Seidel迭代法关于任意初始向量均收敛。其相应的分量迭代格式为:

$\begin{array}{l}\overrightarrow{r}{}_i{}^{\left(k+1\right)}=\frac{1}{{\rm K}{}_{ii}}\left[\overrightarrow{F}{}_i-{\displaystyle \sum _{j=1}^{i-1}{\rm K}}{}_{ij}\overrightarrow{r}{}_j{}^{\left(k+1\right)}-{\displaystyle \sum _{j=i+1}^m{\rm K}}{}_{ij}\overrightarrow{r}{}_j{}^{(k)}\right]\\ (i=1,2,\cdots ,m;k=0,1,2,\cdots )(4)\end{array}$

当边界变形或运动到n+1时刻, 利用 (4) 式经过数次迭代, 便可以达到需要的精度, 得到n+1时刻的网格点的位置。

2 弹簧倔强系数的选择

2.1 标准方法

标准的弹簧近似方法大都选取网格边的边长为参考量来确定弹簧的倔强系数, 大量数值试验表明, 对于二维或三维网格, 弹簧倔强系数按 (3) 式取, 能较好地保持网格的疏密特征。

${\rm K}{}_{\text{s}\text{p}\text{r}\text{i}\text{n}\text{g}ij}=\frac{1}{l{}_{ij}^2}(5)$

(3) 式中lij为网格边的边长。于是当lij→0时, Kspringij→∞。这样便很好地保证了任意的网格节点i, j不会碰撞。

2.2 方法改进

标准的弹簧近似方法仅仅考虑了弹簧在直线方向的伸缩作用, 在二维或三维时, 并不能保证避免网格边的互相交叉。而且, 在用标准的弹簧近似方法处理相对大一些的边界运动和变形时, 就会产生体积为负的“无效”网格单元。为了防止产生这种“无效”网格单元, 就必须在弹簧倔强系数中引入单元几何形状的变量。因此, 在实际计算中对 (5) 式稍作改变, 成为 (6) 式。

${\rm K}{}_{\text{s}\text{p}\text{r}\text{i}\text{n}\text{g}ij}=\frac{1}{V{}_{\min }l{}_{ij}^2}(6)$

(6) 式中Vmin为拥有网格边ij的所有单元中, 体积最小的单元。当Vmin→0时, Kspringij→∞。这样就保证了与网格边ij相对的所有单元中体积最小的单元不会成为体积为负的“无效”网格单元。

为了防止挤压, 已有的四面体单元及四面体中的三角形面元, 必须考虑挤压弹簧倔强系数。挤压弹簧倔强系数计算的具体公式为

${\rm K}{}_{\text{s}\text{q}\text{u}\text{a}\text{s}\text{h}ij}={\displaystyle \sum \frac{1}{\sin {}^2\theta }}(7)$

(7) 式中θ为网格边ij所对应的三角形内角, 以及网格边ij所对应的四面体单元所夹内角, 如图1所示。当θ→0或θ→π时, Ksquashij→∞。这样就避免了挤压已有四面体及三角形面元。

此外, 一些文献把采用弹簧近似方法调整网格点位置的方法, 与椭圆型结构网格生成方法比较。认为, 此种方法也具有椭圆型方法的性质, 即局部的扰动只在局部产生影响。计算也证明了这一点, 最早出现网格变形失败的单元, 就是在运动边界附近。为减小这一性质的影响就必须使弹簧的倔强系数沿内边界到外边界的径向方向合理分布。因此, 我们引进弹簧倔强系数的径向分布系数Kd, 称之为边界修正。

经过上述的改进和修正, 最后得到总的弹簧系数为:

${\rm K}{}_{ij}={\rm K}{}_d({\rm K}{}_{\text{s}\text{p}\text{r}\text{i}\text{n}\text{g}ij}+{\rm K}{}_{\text{s}\text{q}\text{u}\text{a}\text{s}\text{h}ij})(8)$

(8) 式中Kd=1/d, d为计算区域内的点到物面的最小距离。改进后的方法与标准方法的对比结果如图2所示, 图2 (a) 是标准方法的结果, NACA0012翼型绕1/4弦点偏转35度, 尾部已经出现失败。图2 (b) 是对同一问题使用改进后的方法, 此时翼型偏转到50度, 网格质量仍然较好。

3 动网格方法的应用

3.1 控制方程

ALE (Arbitrary Lagrangian-Eulerian ) 有限体积描述下的三维可压缩Euler方程可表示为如下的积分形式:

$\frac{\partial }{\partial t}{\displaystyle \underset{\Omega }{\iiint }\overrightarrow{Q}}\text{d}V+{\displaystyle {\iint }_{\partial \Omega }\overrightarrow{F}}(\overrightarrow{Q})\cdot \overrightarrow{n}\text{d}S=0$

Ω是控制体, ∂Ω是控制体边界, n为控制体边界外法向单位向量, 守恒项和对流项分别为

$\begin{array}{l}\overrightarrow{Q}=[\rho \rho u\rho v\rho we]{}^{\text{Τ}}‚\\ \overrightarrow{F}(\overrightarrow{Q})\cdot \overrightarrow{n}=(U\cdot n)[\rho \rho u\rho v\rho w\rho e+p]{}^{\text{Τ}}+\\ p[0n{}_{x}n{}_{y}n{}_{z}a{}_t]{}^{\text{Τ}}。\end{array}$

其中$\overrightarrow{U}$为流体相对于网格的速度, at为网格运动的法向速度, nx, ny, nz是n的3个分量。

xt, yt, zt分别是网格运动速度的3个分量。空间离散采用格心格式的有限体积法, 时间离散采用Jameson[4]提出的双时间方法, 几何守恒律及刚体动力学模型可参考文献[5,6]。

3.2 算例

3.2.1 简谐振动NACA64A010的跨音速绕流

所计算的机翼剖面形状为NACA64A010, 转动轴定在50%的弦长位置处, 马赫数Ma=0.796, 瞬间迎角α (t) 的运动方程为:α (t) =α0+αmsin (wt) , 式中平均攻角α0=0.0°, 最大振幅攻角αm=1.01°, 减缩角频率K=0.202。为了得到稳定的周期解, 计算了6个周期。

3.2.2 M6机翼变后掠的跨音速绕流

采用的模型是一种以M6机翼为原型的理想的不产生扭转变形可任意拉伸的柔性机翼, 其运动方式是攻角和翼根固定, 翼稍截面向后平移, 机翼自由拉伸, 其展长与弦长不变, 机翼面积也不变。选机翼前缘后掠角为起始角, 在此基础上向后偏转20度, 马赫数Ma=0.80, 总时间T=2.0秒, 攻角a=3.00°。

3.3 结果分析

图3, 图4分别给出了NACA64A010瞬间总升力系数CL和瞬间总力矩系数Cm (力矩作用点为50%弦长处) 随瞬间迎角变化的曲线, 图5给出了在第6周期内8个均分点处的瞬间压力系数分布。图6给出机翼起始和后掠20度后的状态, 图7给出了M6机翼总升力系数随后掠角变化曲线, 图8给出了总力矩系数随后掠角变化曲线。

计算的结果与实验结果吻合良好, 所发展的按弹簧近似的非结构动网格技术与有限体积格式流场解算器相结合, 是一种求包含象机翼变后掠这样的运动边界的三维非定常复杂流场的有效方法。

4 小 结

二维机翼谐和振动结果表明, 本文所发展的弹簧近似方法, 在不降低计算效率的同时, 可以大大提高网格的变形能力和网格质量, 很好地解决了非结构动网格技术中关键的网格变形问题。研究结果表明, 本文所发展的非结构动网格技术及适用于动网格的非定常有限体积解算器具有很高的精度和很强的适用性, 适用于模拟流体与结构耦合等复杂非定常流动问题。

参考文献

[1] Batina J T.Implicit flux-split Euler schemes for unsteady aerodynam-ic analysis involving unstructured dynamic meshes.N91-10918, 1990

[2] Teaduyar T E, Behr M, Liou J.A new strategy for finite elementcomputations involving moving boundaries and interfaces—The defor-ming-spatial-domain/space-time procedure:I.The concept and thepreliminary numerical tests.Computer Methods in Applied Mechacicsand Engineering, 1992;94:339—351

[3] Blom F J.Considerations on the spring analogy.Journal of Aircraft, 2000;32:647—668

[4] Jameson A.Time dependent calculation using multigrid with applica-tions to unsteady flows past airfoils and wings.AIAA paper 1991;91—1596

[5] Hwang C J, Yang S Y.Locally implicit total variation diminishingschemes on mixed quadrilateral-triangular meshes.AIAA Journal, 1993;31 (11) :2008—2015

篇7：多重网格法在非结构网格中的应用

图像去噪是数字图像处理领域中最基础也是最重要的研究课题之一, 对其研究有着很重要的理论意义和实际意义。Rudin和Osher[1]提出的整体变分图像复原模型 (ROF模型) 能在去除噪声的同时保持图像的边缘和细节, 是图像复原中的经典去噪模型, 但ROF模型会在较平坦区域内产生“阶梯效应”。为了克服这一缺点, 改善模型的去噪效果, 陈利霞等人[2]对ROF模型进行了改进, 提出了加权整体变分模型, 并给出了模型的数值解法, 即梯度下降算法。无论是ROF模型还是加权变分模型中, 都含有L1项, 使得数值求解并不容易。为了达到收敛, 梯度下降法需要引进正则化参数, 参数选取越大收敛越快, 但数值解偏离真实解越远。为了解决这一问题和寻求加权整体变分图像去噪模型的快速有效算法, 本文考虑结合对偶方法和多重网格方法对加权整体变分模型进行数值求解。对偶问题与原问题是等价的, 但对偶方法中不需要引进正则化参数, 给数值求解带来了方便[3]。多重网格算法也是求解偏微分方程的最快的数值方法[4]。

1 问题阐述

本文考虑下述加权整体变分 (weight total variation, 简记为WTV) 模型[2]:

其中Gδ*u0表示采用高斯滤波器对带噪图像进行预处理, g (x) =1/ (1+ηx2) , η为大于0的常数。

令α=1/λ, 模型 (1) 可转化为如下约束问题:

运用增广拉格朗日方法可将上述约束问题转化为极大极小值问题:

其中γ为罚参数, p=[p1, p2]为拉格朗日乘子, 利用最优化理论可得:

解式 (2) 可得u=u0-αdiv (p) 相应的对偶方程为

用半隐梯度下降法[5]解式 (3) , 选择τ>0, p0=0, 其迭代 (Chambolle对偶迭代) 格式为:

2 加权整体变分模型的多重网格算法

由于加权整体变分模型的对偶方程 (3) 的不可微性, 用Chambolle对偶算法 (4) 来求解该方程收敛会很慢.尤其是在处理比较大的图像的时候, 图像复原过程耗时非常大。本节提出用多重网格算法[6]求解对偶方程 (3) .

将对偶方程 (3) 的离散格式记成如下形式:

其中

表示最细层的离散解, nL=2L+1-1。NL是非线性算子。由于p1, p2具有零Dirichlet边界条件, 我们假定带噪图像u0∈RnL×nL被嵌入在更大的具有零边界的区域。在第k层, uk有nk×nk个像素值, 其中nk=2k+1-1。在最粗层n1=22-1=3。

我们采用多重网格算法来解非线性方程 (5) 。多重网格算法从最细层L层开始。Rkk-1为k层到k-1层的限制, pkk-1为k-1层到k层的延拓, v1, v2分别为每层的前、后光滑迭代步数。

算法1 (基于WTV模型的多重网格算法) 描述如下:

(1) 参数初始化:取PL0=0, v1=v2=12, μ=1/1502, gL=0, n:=0;

(2) 在层执行一次多网格迭代:

其中PLn是方程 (5) 的解的当前逼近值。从k层到k-1层的一步二重网格迭代Pk=NMGM (Pk, gk, k) 的操作如下:

如果k=1, 精确求解方程N1P1=g1;

否则, 在k层操作:

(3) 计算un+1:

(4) 若满足迭代终止条件, 即:un+1-un燮10-4, 则算法终止;否则, 令n:=n+1, 转向步 (2) 。

这里光滑算子使用迭代 (4) , 但实验发现不能快速收敛。我们将迭代 (4) 修改成如下形式:

其中β为大于0的常数。在最粗层精确求解时, 用迭代格式 (6) 迭代200次。

3 算法的改进

本节首先对算法中的光滑迭代 (6) 进行局部傅立叶分析[6]。令

考虑一般的Fourier分量:

这里。我们把e (k) 1, e (k) 2用含Fourier分量的形式表示为:k

接下来计算光滑迭代 (6) 的维数为2×2的扩大矩阵, 记为Aα, β, 它满足如下方程:

由光滑迭代 (6) 可得到:

这里在局部情况下看做是常数[6]。将上式用有限差分方法离散得:

对上式应用Fourier分析理论, 并取h=1, n=16的网格可得:

这样就可得到扩大矩阵的元素如下:

为了计算Aα, β在高频范围θα, θβ∈[-π, π][-π/2, π/2]的最大奇异值, 我们给定C的值。当τ=1/8时, 结果如表1所示。从表中可以看出, 当C变小时, 光滑速度变慢, 尤其是当C=0.001时, 光滑速度很慢, 几乎不能起到光滑的作用。

我们考虑将光滑迭代 (6) 修改为如下形式:

这里取s为常数。实验发现改进的光滑迭代格式 (7) 能够将多重网格方法提速一半。对改进的光滑迭代格式 (7) 进行局部Fourier分析得到扩大矩阵的元素如下:

计算当C, s取不同值时, Aα, β在高频范围内的最大奇异值, 结果如表3.2所示。从表中可以看出, 当s=-4时效果最好, 而且改进后的Chambolle光滑迭代速度较改进前提高了很多。我们将改进的WTV模型的多重网格算法记作算法2。算法2的光滑迭代和最粗层的迭代采用迭代格式 (7) 。

4 数值实验

数值实验中, 我们比较改进的基于WTV模型的多重网格算法 (NMGM) 和基于WTV模型的Chambolle对偶算法 (CDA) 的去噪的效果和收敛速度。为了定义量化标准用来衡量图像的去噪质量, 我们引入目前常用的图像复原质量评价标准:

峰值信噪比PSNR:

其中, u和u0分别表示复原图像和原始图像, M和N表示图像的长和宽的大小。

实验1:我们用改进的基于WTV模型的多重网格算法 (NMGM) 和基于WTV模型的Chambolle对偶算法 (CDA) 对大小为255255的lena图像进行实验, 添加方差为20的高斯白噪声。此数值实验的参数选取和实验结果如表3所示, 复原图像如图1所示。

实验2:我们用改进的基于WTV模型的多重网格算法 (NMGM) 和基于WTV模型的Chambolle对偶算法 (CDA) 对大小为511511的cameraman图像进行实验, 添加方差为10的高斯白噪声, 此数值实验的参数选取和实验结果如表3所示。

从表3和图1可以看出, 实验用基于WTV模型的多重网格算法和Chambolle对偶算法对不同的图像去噪, 得到的PSNR几乎相同, 复原效果相同, 但WTV模型的多重网格算法所用的时间明显比Chambolle对偶算法少。实验1中WTV模型的多重网格算法所用的时间接近Chambolle对偶算法的1/4, 实验2中接近1/5。

5 结语

本文利用多重网格方法对加权整体变分图像去噪模型进行求解, 通过局部傅里叶分析对算法进行改进。实验结果表明:在去噪效果相同的情况下, 改进的多重网格算法的收敛速度明显优于Chambolle对偶算法。

参考文献

[1]Rudin L I, Osher S, Fatemi E.Nonlinear total variation based noise removal algorithms[J].Physica D:Nonlinear phenomena, 1992, 60 (1/2/3/4) :259-268.

[2]陈利霞, 冯象初, 王卫卫, 宋国乡.加权变分的图像去噪算法[J].系统工程与电子技术, 2010, 32 (2) :392-395.

[3]Chan T F, Chen K, Carter J L.Iterative methods for solving the dual formulation arising from image restoration[J].Electronic Transactions on Numerical Analysis, 2007, 26 (2) :299-311.

[4]黄朝晖, 常谦顺.多重网格插值算子的改进算法[J].应用数学学报.2003, 26 (03) :443-450.

[5]Chambolle A.An algorithm for total variation minimization and applications[J].Journal of Mathematical Imaging and Vision, 2004, 20:89-97.

篇8：多重网格法在非结构网格中的应用

关键词：制造网格,共享,应用服务供应商,资源发现

0 引言

制造网格[1]是网格计算在制造领域内的拓展应用,其本质在于实现资源的全面共享,因此MG引起了许多研究者的关注,他们对MG进行了不同角度的研究。如王爱民等[2]研究了面向制造网格的应用平台及虚拟企业建模,刘丽兰等[3]研究了制造网格中服务质量(Quality of Service,QoS)管理系统,等等。尽管制造网格的研究从不同角度取得了一定的成果,但其共享的本质,即关于如何将企业多余的资源共享出去从而实现企业内外资源的无缝协同却仍未得到彻底解决,它涉及到许多诸如资源的发现、协商、交易、安全等相关问题。总体上制造网格的共享问题可归纳为两方面:一是资源提供,即企业作为整体如何向外发布其多余的资源,从而为协同制造提供物质基础;二是资源应用,即用户如何应用分布在网格上的资源。资源发现贯穿其中,因而成为问题的关键。本文主要对资源的应用进行研究。目前,应用服务提供商(Application Service Provider,ASP)作为一种第三方服务模式,为解决其应用问题提供了途径。但在网格环境下,资源以Web服务形式存在,传统ASP集中提供软硬件资源的服务模式已不再适应网格的应用,尤其是不易实现对异地分散资源的广域发现。许多文献对网格环境下ASP模式如何实现资源发现做了许多有益的研究。文献[4]以ASP平台为网格服务注册中心,ASP节点间采用线性连接,不利于资源的广域发现;文献[5]采用了总控ASP节点与分布式ASP节点相结合的层次结构,实现了多对多的服务模式,各ASP节点存储相邻ASP节点的路由信息,能够实现广域资源的发现,但它仍是通过ASP平台上注册的资源详细信息的细粒度服务信息对资源实现匹配查找,而未从企业服务整体考虑抽象出功能概要等粗粒度信息,从而不能实现快速的分级检索;另外,它需通过ASP总控中心访问网格,查询方法不灵活,而且易形成性能瓶颈。

针对资源应用问题及其现有研究的不足,研究资源广域发现的实现,在文献[5]研究的基础上,提出了多层树型覆盖网络结构的ASP分布,并在该结构基础上对ASP平台进行了功能分析及其设计。

1 多层树型覆盖网络结构的ASP分布

在制造网格环境下,利用ASP向用户提供服务。为了实现多点访问、分层检索以及快速的网格资源发现,对ASP模式的网络结构,在文献[5]基础上加以改进,首次提出构建多层树型网络结构的ASP服务体系,各节点在自治的基础上互相支持。

为了描述方便,首先对该结构中的有关概念进行如下定义:

定义1企业资源节点,是资源所在的网络节点。资源包括设备、软件、智力、人力等。它可以是提供单个资源能力的节点,也可以是提供团体资源能力的节点。

定义2企业信息节点,是企业资源的服务信息管理服务器,它基于通用企业服务中心并配置了其它信息服务相关功能,如资源计费、信息路由等,从而实现资源的注册、发现、交易等信息服务。

定义3 ASP节点,即ASP平台,能实现企业信息节点的粗粒度服务信息的注册,诸如企业信息节点的功能特征、地理范围等概要信息。同时,它能实现用户登录、资源的查询、路由、计费、安全、故障等机制,为用户提供资源共享的平台,并确保网络结构上的信息畅通。

1.1 多层树型覆盖网络结构的ASP分布

ASP分布的网络拓扑结构——多层树型覆盖网络结构如图1所示。

该网络结构以企业资源节点为叶子节点,以企业信息节点为企业资源的信息汇集中心,即注册中心,再向上是更高层的信息注册中心,即其父节点。依此类推,形成多层树型覆盖网络结构的ASP服务体系架构,从而为制造网格的资源发现提供基础。

1.2 多层树型覆盖网络ASP结构特征

1)企业信息节点,它也具有明显的ASP节点的特征。因它负责存储企业资源的注册信息,并对这些资源的使用进行监控并计费管理,同时,它接受外部资源的请求,向外提供资源共享协同服务。

2)除企业信息节点外的各级ASP信息节点不直接注册网格服务信息,仅注册其下级信息节点的粗粒度信息如其功能特征、地理范围等概要信息,为资源的搜索提供基础。

3)各级ASP节点具有严格的层次关系,上级为下级提供信息监控管理服务,而下级为上级节点提供注册信息。下级节点逻辑上是上级节点的子节点,反之,则是下级节点的父节点。

4)各ASP节点即可以下级ASP节点为子节点,也可以企业节点为子节点。

5)各级ASP节点的聚集,即远近距离关系有一定的规律,即满足功能相近、地地理位置相邻的原则,这样构建的ASP制造网格资源发现体系结构能保证搜索到的资源符合制造业的特征,即尽量减少物流运输成本、节省离线作业时间等。

在以上特征中,其中企业信息节点与其它ASP节点相似特征很重要,因为企业信息节点的建设也可由第三方服务商如ASP供应商或其它开发商筹建,两者的相似性致使软件开发能发挥代码重用的功效,从而节省开发成本,降低开发难度,而且由于开发的一致与统一性,方便日常维护。

基于该多层树型覆盖网络ASP结构,用户能通过企业信息节点和企业外的各ASP节点实现资源的发现查找,从而实现广域多点发现。另外,各ASP节点相互支持,但又彼此保持高度的自治性。同时,通过ASP节点上的企业级或ASP级别的粗粒度信息的检索定位,实现了基于层次的分级检索,有利于快速定位。

2 资源发现并行搜索过程

在基于上述资源发现结构支持下,任何信息节点都可以处理节点的资源查找请求。当用户或应用(这里统称为用户,以下同)查找所需资源时,其资源的查找过程如下:

1)首先用户根据注册帐号登陆相应的ASP节点。用户分两类:一类是企业内用户,选择从企业内部登陆,其资源搜索路径如需向企业外部搜索则需经该企业信息节点;另一类是企业外用户,选择从外部注册ASP信息节点登陆。

2)企业内用户通过企业内部用户帐号从网格门户的企业节点入口登陆到网格系统,向企业提出资源请求,然后进入下一步,如果是企业外用户则转6。

3)在企业节点历史库或注册资源库中,资源请求进行资源匹配查找(在通用企业服务控制中心实现),如果找到符合要求且数量足够的资源则一次性返回资源信息,然后进入下一步;否则企业节点将查找请求发送到上级ASP节点,然后转5。

4)用户端根据一系列资源选择策略自动进行择优筛选,如考虑资源的Qo S,资源的信息路径长度,即能选择短链的则不选择长链等,进行资源链组合,得到最优资源链,然后用户发起协商请求,与服务提供者协商并签署协商协议,确定的资源信息路径在信息返回中被各个路径节点分别根据自身位置记录其相应内容,为计费提供依据,也为日后查询提供缓存数据(资源的信息被分类地记录在各确定的信息路径节点上,显然高层节点将是对其各自节点的历史汇总),然后转7;若未达成协议,为避免重复搜索,根据所确资源种类、数量继续从先前的搜索位置的下一位置开始搜索,如能找到足够满足数量的资源则不需向上层搜索而返回,否则,将搜索层次加1,进入下一步。

5)ASP节点根据其注册资源的功能特征选择其搜索路径,如该ASP节点没有符合的功能特征,则继续向上搜索,直到找到符合需求功能特征的层次为止;若有则先在资源历史库中搜,并判断资源数量是否符合数量要求,够量则返回,反之则并发的将请求发往该层符合要求的各子信息节点,收到消息的节点又分别以并行方式向各自的子节点发送资源请求,直到企业信息节点,然后转3。

6)企业外用户通过注册帐号从网格门户选择注册ASP节点入口登陆到网格系统,向该节点发出资源请求,然后转5。

7)资源发现停止。

该网络结构上的资源查找是整个资源发现的关键,分为路由和目标节点的匹配查找两个阶段。不失一般性,可将树型结构看作是均匀组织的。从上述过程可知,基于该结构资源搜索采用并行搜索算法,即各子树并行搜索。

3 ASP平台的设计

3.1 ASP平台功能需求

根据上述ASP分布结构及其资源查找过程,ASP平台的所需主要功能要求如下:

接受用户的注册与登陆;接受用户的查询请求,并对请求进行解析,分解成面向服务的查询请求;对查询请求,首先在历史库中进行服务匹配,未找到所需的服务则在其注册库中进行分类查找,根据其功能、地理范围类别找到下级资源路径;根据路径,系统对其进行路由,将请求发向下级节点;在信息更新变动中,ASP节点需及时维护对资源的注册变化信息,确保信息的一致性、实时性;对ASP节点所在的上下文环境中,需对其子节点进行监控,对父节点进行及时上报,确保其信息畅通;对找到的资源进行最佳资源链配置以及ASP平台的基本管理功能需求如计费管理、用户管理、日志管理、安全管理,等等。

3.2 ASP平台软件设计

根据上述ASP平台的功能需求分析,参考文献[6]的ASP平台软件结构图,得到本文研究的ASP平台软件体系结构,如图2所示。各层的功能如下:

数据层:负责综合管理ASP平台中的信息资源,提供信息的存储与共享服务,其中用户信息数据库是存储用户的注册、访问权限等信息,历史数据库存储该ASP节点的以往资源查询信息等。

平台管理层:主要为平台提供服务环境,包括用户管理、安全管理、日志管理、运行监控、计费管理、等常见功能模块。其中计费管理主要包括对本节点所辖资源的使用计费以及充当资源信息路径的管理费;运行监控是该节点对其父子节点运行情况的实时监控;用户管理是对该节点的注册用户的信息管理,是用户登录并获取资源的重要信息依据。

使能资源层:该层是构建面向用户的应用系统的资源集合,它为用户的应用提供资源构建基础,最终向用户呈递基于制造网格系统的优化的资源链配置,致使用户通过ASP平台能用到满意的资源,得到满意的结果。其中UDDI注册中心依据分类的原则对来自ASP节点内外的资源进行注册,并提供基于分类的资源查询;资源注册与查询提供基于UDDI注册中心的信息注册与查询功能;信息更新与维护提供对ASP节点内外的资源变化的更新与维护。上述三部分类似于制造网格资源发现结构的企业信息节点的信息服务中心的相应部分,这样可减轻开发的负担,达到代码重用的目的。信息路由是将本ASP节点的信息查询请求按照所选路径向外部ASP节点发布,或是将信息反馈按照指定路径向下一节点进行转发等,由此使ASP节点在资源发现结构中起到中间桥梁的作用,从而帮助实现资源的广域发现。资源链优化配置是ASP节点对查询到的资源进行资源链优化配置,为用户选择满意的资源,从而为用户提供满意的服务并使他们得到满意的结果。

应用系统层:由于不同用户具有不同的需求,需要根据他们的要求基于资源使能层构建适合其需要的面向服务的应用系统。由于ASP平台的广域发现性,因此能构建不同应用的基于服务的系统来。

用户界面层:该层是直接面向用户的。用户通过它登陆系统,向系统提供资源查询请求,并通过它得到资源反馈结果,还能通过它了解到任务的执行情况等。该层建设的好坏,直接影响到系统的使用,是该平台建设乃至整个制造网格系统建设的关键一环。

该软件体系结构的特色如下:

由于本文的ASP平台是面向网格服务,它不是面向实际的资源,因此该ASP平台的内部功能元素是面向网格服务的,而没有针对实际物理资源的管理功能,致使ASP平台的通用性较强,可以接受来自不同用户的不同制造请求;

该ASP平台既能实现对异地资源的调用,又能实现对本地资源的调用,因而致使其实用性较强,也便于对平台的功能扩充;

另外,由于该ASP平台有别于传统的集中式管理的模式,它既能充当集中管理资源的平台(对本地资源)又能充当信息路由的中间节点,因此是一种集中与分布式结合的管理方式。

这里重点对上述ASP平台软件体系结构的计费管理及信息路由进行论述。

计费管理:结合现行ASP模式收费机制特征及该模型的特殊性,各ASP节点收取费用可通过以下渠道实现:1)各ASP节点向其直接下级节点收取注册年费,作为为下级节点提供信息维护、信息监控等的日常维护开支费用;2)跨域收费。根据地理范围判断其跨域情况,并按跨域收费标准收取费用;3)信息服务费,每次信息查找,资源信息路径是确定的,当成交后,该条信息路径中的各信息节点记录本次交易路径信息作为计费依据,作为从其下属节点收取信息服务费用的依据。而企业方根据客户的信息路径及其制造时间等为依据来决定收取的资源价格。为防止客户可能会选择靠近信息节点的位置登陆,从而降低其交易费用,各节点都有注册帐号使用限制,即何处注册何处登陆,并且系统限制单人单账号,杜绝个人以不同身份不同点登陆,另外信息反馈只能反馈到登陆点。为防止客户欺骗行为,企业方必须据客户信誉度达成协议,若成交,各信息节点向其下级节点收取每次交易信息服务费用。

信息路由:从一ASP节点发出资源请求到发现资源的算法,由资源发现发起算法、资源路由算法、企业内资源查找算法等组成,信息路由算法伪码描述如下:

4 实现

基于该模型,用GT4中间件构建了制造网格平台,并选用开源工具JUDDI来实现扩展的企业UDDI注册中心;各ASP信息节点初步采用目录结构存储方式实现了简单的概要信息记录;基于文中的信息路由和其它相关的资源发现算法,并基于一简单的ASP客户端界面实现了对简单资源的访问。图3是企业信息节点的原型系统实现框架。

5 结论

本文从资源发现的角度,针对ASP集中分布式结构的不足,提出了多层树型覆盖网络结构,并设计了ASP平台的软件体系结构。以上研究对制造网格的应用具有积极的意义。下一步,我们将继续完善ASP平台的构建工作,进而对资源的调度、安全等方面开展深入的研究。

参考文献

[1]范玉顺,刘飞,祁国宁.网络化制造系统及其应用实践[M].北京:机械工业出版社,2003.

[2]王爱民,范莉娅,肖田元,等.面向制造网格的应用平台及虚拟企业建模研究[J].机械工程学报,2005,41(2):176-181.

[3]刘丽兰,俞涛,施战备.制造网格中服务质量管理系统的研究[J].计算机集成制造系统,2005,11(2):284-288.

[4]牛秦洲,叶恒舟.基于交互式网格的ASP体系结构[J].电信科学,2005,4:72-73.

[5]王国庆,王刚,吕民,高国安.基于网格的应用服务提供商平台制造资源共享方法研[J].计算机集成制造系统,2007,13(2):350-355.