徐兵兵团队《IJNME》新作：虚拟元法求解压电接触问题（附开源程序）
《International Journal for Numerical Methods in Engineering》最新论文：
Virtual Element Method for Piezoelasticity
近期，徐兵兵博士后作为通讯作者，和国内学者合作的最新虚拟元法(Virtual Element Method,VEM)的论文正式上线。该文章使用VEM求解了经典的多物理场耦合的压电问题，并考虑了接触非线性，并附有开源代码。
论文原文见：https://doi.org/10.1002/nme.70191
开源代码见：https://www.vemhub.com/1058.html
摘要
本文提出了一种用于模拟二维和三维压电问题的虚拟单元法（VEM）。压电材料具有很强的多物理场耦合特性，能够将机械能转化为电能。由于需要在统一的框架内处理力学、电学和接触约束，因此精确的数值模拟极具挑战性。VEM 是一种新的数值方法，能够处理一般的多边形和多面体网格，使其适用于模拟压电问题以及涉及接触的情况。其核心思想是从本构方程中分解总能量密度，从而推导出耦合机械-电问题的双线性公式。本文还讨论了耦合问题的稳定项。在 VEM 框架中，二维情况下的接触采用罚函数法进行处理。
基本理论
变分形式
首先引入两个关于位移和电场的函数空间
对于压电问题，其总势能为
紧接着，考虑其变分 ，可得到
最终，我们得到以下关于压电问题的耦合方程
其中
虚拟元法
在该工作中，考虑了二维及三维问题。在虚拟元法中，其核心思想为投影算的的定义及计算。 考虑任意多边形单元，位移场变量 （或者电场 ）可以通过投影分为以下两部分
其中的 投影算子可通过以下正交形式计算
对于二维单元，我们采用以下自由度的定义
在单元角点的节点值 位于边界 上的 个Gauss-Lobatto点的值 单元内部矩量
通过以上自由度，我们可以求解投影算子 . 对于三维问题，该工作仅考虑一阶格式。
离散格式
通过计算VEM中的投影算子，则对于压电问题的耦合变分方程，可得到其离散形式为
其中 为稳定项
考虑投影算子的矩阵形式，其单元刚度矩阵为
对于该耦合场问题，稳定项参数可选择为
压电接触问题
对于许多压电问题，其多需要考虑外接接触带来的驱动。对于经典的接触问题，可以采用拉格朗日乘子法或者罚函数法处理该约束优化问题。在VEM中，我们采用了相同的技术。但是由于VEM可以采用多边形单元，因此可以使用点-点接触进行接触面的离散。 考虑罚函数法，其接触势能的变分为
则接触刚度矩阵为
数值算例
使用多边形网格的二维压电接触问题
三维压电问题
开源代码
本工作同时开源了VEM求解二维压电问题的MATLAB代码，可以通过VEMhub的官网或者github网站链接下载。具体的下载链接为
