RKDG方法改进的K均值聚类坏单元指示子

赵郑豪, 王之欢, 朱洪强

(南京邮电大学理学院, 南京 210023)

非线性双曲守恒律方程的数值解通常会出现间断,导致数值计算中产生振荡, 造成不稳定性, 给数值模拟带来很大的困难, 而RKDG方法(Runge-Kutta discontinuous Galerkin methods)[1]是求解这类方程的主流方法之一, 其主要过程为: 1) 通过设计坏单元指示子快速定位间断位置, 并标记出附近的坏单元; 2) 重构坏单元上的数值解控制数值振荡．RKDG方法的另一个特点是近似解不要求单元间的连续性, 故该方法能在各类网格上使用, 包括带有悬点的非协调网格．目前已有许多关于坏单元指示子的研究工作, 在较早的工作中, Qiu等[2]数值比较了7个常用的坏单元指示子, 但并没有得到一个对多数问题表现一致最优的坏单元指示子; Dumbser等[3]对2014年之前的坏单元指示子工作做了比较详尽的回顾．近几年, Vuik等[4]提出一种根据不同坏单元指示变量能自动选择参数的坏单元指示子; Gao等[5]在此基础上进行了拓展优化; Fu等[6]提出一个既简单又不含敏感参数的坏单元指示子; 笔者[7]则将此工作推广到了h自适应网格上．近些年, 机器学习在这一领域也得到了广泛应用, 如基于神经网络设计的坏单元指示子[8-11]等．

间断区域和连续区域内单元的指示变量值往往相差较大, 笔者[12]基于这个观察设计了一个基于K均值聚类的坏单元指示子, 该工作通过将所有单元划分为若干个局部小模板, 让每个模板中至少包含一个好单元, 再利用统计方法中常用的K均值聚类方法将模板中的单元最终分成两类, 然后使用一个新设计的模板分类准则来判断模板中是否存在坏单元, 从而确定每一个单元是否为坏单元．这个K均值聚类坏单元指示子对局部收集来的指示变量值运用数据分析方法做坏单元判断, 避开了对网格剖分的依赖, 可以用于各种类型的网格, 且在数值试验中能准确判断间断位置．然而,此方法含有2个实数参数, 虽然它们不依赖测试问题, 但参数的确定需要大量的数值测试, 影响了该方法的实际可用性．为了解决这个问题, 本文拟基于K均值聚类坏单元指示子[12]的框架, 在聚类之前首先对指示子变量值进行归一化来消除不同算例中指示变量取值范围不同的影响,然后利用极差判断模板中的单元是否存在坏单元．这种做法不仅简单、高效,而且仅含1个参数, 能极大地降低确定参数的复杂度, 有利于该方法在实际工作中的应用．

受文献[12]坏单元指示子框架的启发, 本文将坏单元指示子仍然分为两部分: 1) 为了让探测坏单元在解结构更为简单的局部区域进行, 将计算单元按一定规则划分为一些模板; 2) 对每个模板进行坏单元探测．具体步骤:

i) 对于第(1)部分, 模板划分的方法参见文献[12]．因为全是好单元和全是坏单元模板上指示变量的值在数据表现上不易区分, 故要求每个模板中至少有1个好单元．同时, 为了节省计算时间, 模板大小(即模板中包含的单元数)越小越好, 故一维时模板尺寸取为5, 二维矩形网格时模板尺寸取为4×4[12]．

ii) 对于第(2)部分, 需要对每个单元判断是否为坏单元．这部分又可细分为: 归一化处理、模板筛选和坏单元筛选．由于坏单元指示变量对于不同的算例或不同的DG多项式次数k, 其取值范围各不同, 导致坏单元指示子中的参数可能会依赖于算例或k值, 给数值应用带来麻烦．受文献[13]的启发, 首先对指示变量的值进行归一化处理．假设计算单元总数为n, 给定指示变量I, 所有单元上的指示变量值分别记作I1,I2,…,In．定义Imax=max{I1,I2,…,In},Imin=min{I1,I2,…,In}．对每个单元的指示变量值Ii做线性变换

Ti=(Ii-Imin)/(Imax-Imin)

,

(1)

指示变量须进行归一化处理, 否则对不同情况下参数σ的最优取值明显不同[13], 这给实际应用带来不便．经归一化处理的坏单元指示子可使用多种坏单元指示变量, 本文采用单元边界跳量指示变量[14]．假设目标单元为K0, 且有G个直接相邻单元, 分别是K1,K2,…,KG．若pj(j=0,1,…,G)是单元Kj上的DG近似解, 则单元边界跳量指示变量

注与文献[12]的做法相比, 本文算法中去掉了原来含有2个参数的过滤步骤和模板分类判断步骤, 改为使用归一化来优化数据, 结合只含有1个参数的极差法做模板分类判断, 能简单、高效地筛选出含有坏单元的模板．由于新的坏单元指示子只涉及1个参数, 少量的测试就能确定参数取值, 极大地方便了该法的数值应用．

以欧拉方程组作为模型方程, 采用经典的算例来展示改进的坏单元指示子的效果．若方程在二维情形时的形式为

其中ρ为密度,E为总能量,p为压力,u和v分别为x和y方向上的速度分量, 这些量满足p=(γ-1)[E-0.5ρ(u2+v2)], 式中γ=1.4．改进的坏单元指示子基于密度和能量分别探测坏单元, 一个单元只要有一次被探测为坏单元, 最终就会被标记为坏单元．

例1Lax问题．有初值条件

取网格单元大小为Δx=1/40, 图1(a)和(b)为t=1.3时的密度数值解, 图1(c)和(d)为不同时刻坏单元的位置图(坏单元历史图)．计算结果表明, 改进的坏单元检测方法能准确标记出间断的位置, 数值解没有产生振荡．

图1 Lax问题在t=1.3时的结果Fig.1 Results at t=1.3 for the Lax problem

例2爆炸波问题．初值条件为

图2给出了Δx=1/800的网格在t=0.038时的数值解．从图2可以看到, 改进的坏单元指示子准确探测到了2个大间断, 且数值解无振荡．

图2 爆炸波问题在t=0.038时的结果Fig.2 Results at t=0.038 for the blast-wave problem

例3双马赫反射问题．设计算区域为[0,4]×[0,1], 底边1/6≤x≤4区域处为固壁．当过点(1/6,0)且与固壁成60°的激波以10马赫的速度沿垂直方向撞击固壁时, 设激波前方没有受到扰动的空气密度为1.4, 压力为1; 底侧的固壁区域使用反射边界条件, 其他区域使用精确的波后边界条件; 顶部边界使用10马赫激波运动的精确边界条件, 左右边界依次使用入流边界条件和出流边界条件．图3给出了t=0.2,Δx=Δy=1/240时得到的密度等高线和坏单元．结果显示, 虽然在k=1的情形中间断位置没有被完全探测出来, 但所有的数值解没有振荡．改进的坏单元指示子在确保数值解没有振荡的前提下,标记了尽可能少的坏单元．

图3 双马赫反射问题在t=0.2时的结果Fig.3 Results of the double Mach refection problem at t=0.2

例4前台阶问题．问题起始于一个长度为3, 高度为1的风洞中, 风洞底部有一高为0.2、离风洞左边端口0.6并一直延伸到风洞右端口的台阶; 风洞壁处使用反射边界条件, 左右边界分别使用入流和出流边界条件．图4给出了t=4.0, Δx=Δy=1/160时得到的密度等高线以及坏单元．从图4可以看到, 坏单元探测准确, 数值解并没有振荡．

图4 前台阶问题在t=0.4时的结果Fig.4 Results of the forward-facing step problem at t=0.4

为了验证改进的方法仍然适用于不同的坏单元指示变量, 图5给出这个算例使用Fu-Shu指示变量[14]的计算结果．从图5可以看到, 坏单元检测准确, 数值解没有振荡．

为了说明改进的坏单元指示子在不同粗细网格上的表现, 表1列出了这个算例2个指示变量在两套网格上的坏单元平均百分比．结果显示, 当网格尺寸减半时, 坏单元平均百分比在所有情况下都大约减小到原来的一半, 说明改进的坏单元指示子在不同粗细网格上的表现良好．

表1 前台阶问题坏单元平均百分比

猜你喜欢算例边界条件聚类一类带有Stieltjes积分边界条件的分数阶微分方程边值问题正解应用数学(2020年4期)2020-12-28带有积分边界条件的奇异摄动边值问题的渐近解数学物理学报(2020年5期)2020-11-26基于DBSACN聚类算法的XML文档聚类电子测试(2017年15期)2017-12-18基于高斯混合聚类的阵列干涉SAR三维成像雷达学报(2017年6期)2017-03-26基于振荡能量的低频振荡分析与振荡源定位（二）振荡源定位方法与算例中国学术期刊文摘(2016年2期)2016-02-13互补问题算例分析新乡学院学报(2015年6期)2015-11-06基于CYMDIST的配电网运行优化技术及算例分析电网与清洁能源(2015年2期)2015-02-28带Robin边界条件的2维随机Ginzburg-Landau方程的吸引子四川师范大学学报（自然科学版）(2015年1期)2015-02-28一种层次初始的聚类个数自适应的聚类方法研究电子设计工程(2015年6期)2015-02-27燃煤PM10湍流聚并GDE方程算法及算例分析电力工程技术(2014年5期)2014-03-20

推荐访问:指示 单元 均值