自适应有限元技术:一维电场衰减数值模拟

简述

    在前几篇文章详细介绍了一维有限元的各种实现方法与过程，这里再介绍下自适应有限元技术。

    自适应有限元技术的目的：通过不断修正网格粗细或者基函数阶数，以获得数值模拟的最精确解。其主要包含三个部分：有限元数值模拟、后验误差估计、网格加密策略，这三个部分又对应了不同的解决方案，例如有限元数值模拟分为线性、2阶、高阶等等；后验误差估计分为基于残差、基于恢复梯度等方式；网格加密自然就是不同的网格剖分算法。

    这里通过介绍一维自适应有限元的技术，完整的介绍了自适应有限元技术的整个流程。使用的边值问题：一维电磁波在导体介质上的衰减规律。该问题的有限元实现参考：一维混合高阶有限元详细实现过程。这里不再过多赘述。

1.通过几个例子，了解为什么要使用自适应技术

    使用之前文章一维混合高阶有限元详细实现过程 的代码计算20个网格单元下的不同频率电场衰减变化。具体参数如下：



    当网格均匀且数量不变的情况下，随着频率变化，网格中心的电场误差变化情况如下：


a.100Hz,网格量=20




b.1000Hz,网格量=20



c.10000Hz,网格量=20




  通过上述例子发现，每个频率误差较大的地方均出现在左侧，随着频率增加，误差在逐渐增大。在低于100Hz时候，当前网格还能满足计算精度，误差低于2e-3，但是当频率高于1000Hz的时候，误差已经高达0.04，相对于绝对误差已经是非常大了。




    从中我们发现两点问题：

     I.低频的时候只需要较少的网格就能实现高精度，高频的时候同样的网格不再满足精度。

    II.所有频率在左侧的计算精度明显要低于右侧的计算精度，也就是右侧的网格满足精度要求，左侧的网格不满足精度要求。




    针对问题，把网格粗暴加密是一个直接了当的办法，例如：

a.100Hz,网格量=200



b.10000Hz,网格量=200



    可以发现，在使用200个网格单元后，即使是高频也能保证计算误差小于2e-3以下。

    但是同样会发现三个问题：

        I.低频完全不需要这么多的网格；

        II.网格右侧的计算精度同样不需要如此密的网格。

        III.如果频率继续增加，是否又必须人工增加网格个数。




    因此为了解决上述的问题，自适应有限元技术就孕育而生。

    自适应技术就是根据不同边值问题，计算出后验误差估计，然后将误差大的地方进行加密网格，误差小的地方不加密或者放粗网格，最终获得计算精度解与计算资源共同实现最优的效果。

   自适应有限元技术的流程大致如下：






    在整个流程中，后验误差估计与网格优化是自适应技术的两个重点。这里主要介绍后验误差中恢复梯度的原理以及网格优化方法。




2.后验误差

    摘自一篇论文的描述：




    对于一维1阶有限元而言，同样满足上述两条规律：在节点处计算精度最高，在高斯点（线单元中心点）解的梯度精度最高。




    下图详细示意了第二点规律。在频率1000Hz下，数值计算的电场梯度与理论电场梯度的对比结果。

    明显看出由于采取的是线性1阶基函数，因此在每个单元内的数值梯度是一个常数，在高斯点上的数值梯度的结果是最接近理论梯度的。






    

    恢复梯度的原理就是通过这些高斯点的精确梯度解来恢复其他误差大的地方的精度。

    对于一维1阶有限元而言，在已知高斯点的梯度后，恢复每个节点位置的梯度，最简单的方法就是线性插值，把精确解插值到节点上。









    上图蓝色圈点就表示在节点上恢复后的梯度解，整体上的精度是高于原有梯度解，因此通过对二者的对比，就可以获得整个区域关于电场梯度的误差分布的结果，这个结果就是所谓的后验误差估计。




3.网格加密策略

   在获得后验误差估计以后，就可以根据需求制定精度要求，也就是退出自适应的标准，选择方法也是千奇百怪，本文选择的方法是如下：

 



     

    其次，还需要设置网格加密程度系数，直白的意思就是每次加密多少个网格。系数范围（0,1]，最小表示不加密网格，最大表示全部加密原有网格，这里设置0.5则是加密原有网格数量的一半。    

    

    最后就是如何加密网格，对于1维而言相对简单，直接选择需要加密网格的高斯点，也就是中心作为加密点即可，新的单元则在原单元上一分为二，变成两个单元，例如示意图中加密原有2号单元。






    由于1维模型简单，所以在加入点后，对所有点进行排序，以此重新获得新的网格。

    这种简单的处理方式在2、3维度是非常麻烦的，高维的处理方式一般是保留原有网格的不变的拓扑关系，把新加密点添加在原有网格之后，对网格单元只改变拓扑关系有变化的网格。





4.自适应测试与分析

    在设置好基本参数后，就可以开始进行完整的自适应有限元计算。






a.1000Hz，初始网格=8







b.1000Hz，自适应迭代4次



c.1000Hz，自适应迭代9次



d.1000Hz,自适应误差收敛曲线



    从初始网格、第4次迭代和最后一次的迭代过程中看出，自适应加密的区间的确是误差较大的区域，最终呈现左侧密集于右侧的网格趋势。整体误差低于1.5e-4。自适应收敛曲线在前几次的加密过程中收敛较快，后面几次收敛逐渐放缓。基本上整个自适应流程符合预期。

    下面再看看文章开头的提及的几个频率下的自适应结果：

e.100Hz自适应迭代10次收敛结果



f.10000Hz自适应迭代8次收敛结果



g.100000Hz自适应迭代10次收敛结果



    结果显示自适应结果均满足预期要求，即使频率在1e5Hz，电场在最右侧急剧衰减的情况下，依旧能保证自适应收敛到足够精度。再对比每个频率迭代信息：



    所有频率的最大误差均小于1e-3,这是最开始的200个均匀网格做不到的。

    同时发现频率越高，最后迭代的单元数量反而是更少的，这是由于越是高频，需要加密的地方越集中在左侧，右侧大量的区域不需要加密。例如1e5Hz频率下，在整个0~1000范围内，200~1000的范围几乎没有加密。这也再次说明了自适应技术的必要性。




5.总结

（1）自适应有限元技术的本质目的是为了使用更少的计算资源获取足够的计算精度，这在大规模复杂的三维例子显得尤为重要，好的自适应技术能在有限的计算资源下快速收敛获得精确结果。




（2）万变不离其宗，自适应有限元技术的流程：

        a.给出初始网格、物理参数、自适应迭代加密相关系数；

        b.进行有限元数值模拟；

        c.计算后验误差估计,判断精度是否满足；

        d.根据后验误差设定加密策略，优化当前网格；

        e.重复b~d流程，直到满足精度需求，退出自适应迭代。