第39卷第5期力学与实践2017年10月Trefftz有限元法的研究进展
王克用1)李培超
(上海工程技术大学机械工程学院，上海201620)
摘要Trefftz有限元法(Trefftzfinite element method,TFEM)是一种高效的数值计算方法，兼有传统有限元法和边界元法的诸多优点.基于双独立插值模式，结合杂交泛函和高斯散度定理，推得仅含边界积分的
有限元格式.简述了过去10年间(2007—2016)Trefftz有限元法在单元域内插值函数、源项处理、特殊功能
单元以及非各向同性材料等方面的研究进展，并对未来的发展趋势给出了几点展望.
关键词Trefftz有限元法，域内插值函数，边界积分，特殊功能单元，无源化处理
中图分类号：O343.1文献标识码：A doi：10.6052/1000-0879-17-115
RESEARCH ADV ANCES IN THE TREFFTZ FINITE ELEMENT METHOD
WANG Keyong1)LI Peichao
(School of Mechanical Engineering,Shanghai University of Engineering Science,Shanghai201620,China)
Abstract The Trefftzfinite element method(TFEM)is an efficient numerical approach with many joint advan-tages of the conventinalfinite and boundary element methods.Based on the mutual independent interpolation modes,thefinite element formulation involving the boundary integrations only is derived by incorporating the hybrid functional and the Gaussian divergence theorem.The research advances in the internal interpolation function,the treatment of the source term,the special-purpose element and the nonisotropic material during the past decade(2007—2016)are reviewed and several directions are pointed out for the future development.
Key words Trefftzfinite element method,internal interpolation function,boundary integration,special-purpose element,non-source treatment
引言
Trefftz方法是求解偏微分方程的高效数值方法之一，它是以德国数学家Erich Trefftz(1888—1937)的名字命名的，以纪念他在数学领域的开创性贡献，但在当时并未引起广泛重视.在分析网格畸变对薄板单元的影响时，Jirousek等[1]于1977年首次提出了杂交Trefftz有限元模型.该模型在单元内部和边界上假设两套独立的位移场：单元域内场和网线场.单元域内场可认为是Trefftz有限元法区别于传统有限元法的重要标志，其插值函数(常称为Tr-efftz函数)精确满足问题的控制方程.相邻单元通过网线场在某种杂交意义上关
联成离散的有限元模型.为了推动Trefftz方法在计算力学领域的发展及应用，自1996年开始，每3年召开一次国际会议，报道Trefftz方法的最新进展及研究成果.表1列出了历届Trefftz方法国际会议情况.基于Trefftz方法的数值方法主要包括Trefftz有限元法[2-3]、Trefftz 边界元法[2,4]以及Trefftz无网格法[5-7]，其详细论
2017–04–05收到第1稿，2017–05–03收到修改稿.
1)王克用，博士，副教授，主要研究方向为Trefftz有限元法和多孔介质传热.wang@126
引用格式：王克用，李培超.Trefftz有限元法的研究进展.力学与实践,2017,39(5):433-440
Wang Keyong,Li Peichao.Research advances in the Trefftzfinite element method.Mechanics in Engineering,2017, 39(5):433-440
434力
学与实践2017年第39卷
表1历届Trefftz 工程计算方法国际会议
届次年份
地点内容/主题
truncate的特征11996波兰(克拉科夫)Trefftz 工程计算方法最新进展和展望
21999葡萄牙(辛特拉)Trefftz 工程计算方法32002
英国(埃克塞特)
计算流体动力学中的Trefftz 工程计算方法和有限元法
42005斯洛伐克(日利纳)Trefftz 工程计算方法52008比利时(鲁汶)Trefftz 工程计算方法62011中国(台湾)Trefftz 工程计算方法、基本解工程应用方法
7
2015
中国(杭州)
Trefftz 工程计算方法、基本解工程应用方法
述可参见相关著作[2-3,5-6].
与传统有限元法相比，Trefftz 有限元法的显著优势表现在以下几个方面：
(1)单元刚度方程仅涉及边界积分，这样可对曲线或多边形几何边界进行建模[2-3,8-9]，同时也使得Trefftz 单元对网格畸变不敏感[10-12]，降低了求解维
数，从而减小了计算量.
(2)单元域内插值函数无需精确满足单元间的连续性，而是通过变分泛函使之近似满足[13-16].
(3)通过寻求局部解函数构造域内插值函数，
可捕捉不连续载荷、裂纹、孔洞、夹杂等局部效应问
题[2-3,17].
经过40年的发展，Trefftz 有限元法已经从平面弹性问题、板弯曲问题、位势问题、压电问题等领域拓展到一些新领域，如牛顿流体的流动问题以及软组织[18-19]或水饱和多孔介质问题[20-22].本文仅简述过去十年来(2007—2016)Trefftz 有限元理论的研究进展.
1单元域内插值函数的构造原理
Trefftz 有限元法在单元域内和单元边界上假定
两套独立的场函数插值模式(图1).单元域内场可表达为
u e (x )=
u e (x )+N e (x )c e (x )(在Ωe 内)
(1)
而单元边界网线场为
˜u
e (x )=˜N e (x )d e (x )(在Γe 上)(2)
式中，Ωe 为Γe 围成的单元区域，
u e 为控制方程中
源项诱发的特解，u e 和˜u
e 分别为单元域内和边界上的场变量，N e 为Trefftz 插值函数，是控制方程
的齐次解，˜N e 为网线插值函数，可按传统有限元方法构建，c e 和d e 分别为待定参数和单元节点自由度列阵.m 为截断的完备解或基本解虚源点的个数.
图1基本解Trefftz 单元及其源点
在Trefftz 有限元理论中，基本形成了基于完备解和基本解的两种单元构造方法，所构造的单元可分别称为完备解单元和基本解单元.完备解单元出现得较早，其严格的数学理论是由Herrera 及其合作者完成的[23-26].满足控制方程的完备解系有无穷多项，在构建域内插值函数时，需要适当截断而得到有
限项(亦称Trefftz 项数)[2].为了避免出现零能模式并使单元性能稳定，用来构建N e 的截断完备解个数必须满足不等式m  n d −n r ，这里n d 和n r 分别为单元节点自由度和刚体运动模式的数目.值得注意的是，此不等式是一个必要但非充分的条件.一般地，Trefftz 项数取最小值并不能保证单元刚度矩
第5期王克用等：Trefftz有限元法的研究进展435
阵满秩.在实际应用中，通常需要选取更多的Tre-fftz项数以确保所构造的单元性能稳定，但另一方面项数过多又会导致数值溢出或者病态系数矩阵.因此，针对不同的工程和物理问题，需要数值验证确定合理的Trefftz项数.Choo等[11]利用Mindlin-Reissner厚板问题解构造出9自由度三角形(T32-7)和12自由度四边形(Q32-11)板弯曲单元，当这两种单元退化到薄板极限情形下不会出现自锁现象. Moldovan等[20]基于Biot理论，利用Navier控制方程的自由场解构造域内试函数，从而建立了饱和多孔介质的杂交Trefftz应力元和位移元模型.Wang 等[12]从准调和多项式出发，推得轴对称位势问题的完备解插值函数，并构造了4节点四边形轴对称环状单元.Rezaiee-Pajand等[27]分析了薄板弯曲问题，并构造了两种高阶杂交Trefftz单元：三角形单元(THT-15)和四边形单元(QHT-23).为了更好地兼容，采用3节点Euler-Bernoulli梁的形函数构造网线函数.
由于很难得到一些工程和物理问题的完备解，且截断Trefftz项数时需要格外小心，才能达到预期精度.为了克服这一缺点，基本解Trefftz单元应运而生[28].尽管基本解Trefftz单元采用与传统边界元法相同的基
本解形式，但二者有本质区别.根据互易定理，传统边界元法涉及边界积分方程，在处理奇异或超奇异积分时会遇到困难，而基本解Trefftz单元因不涉及边界积分方程完全不受其限.基于杂交Trefftz有限元思想，Qin[2]和Wang等[28]最早提出了基于基本解的Trefftz有限元模型，并成功应用于二维单层和多层材料的热传导分析中.随后，利用基本解Trefftz有限元法，Wang等[29]分析了正交各向异性平面弹性问题.Cao等[30]利用格林函数形式的基本解作为域内插值函数分析了平面压电问题，并且对应力集中现象提出了一些新见解.Wang 等[13]分析了Poisson-Boltzmann方程和扩散反应方程等两类二维Dirichlet问题.首先在通过每个Pi-card迭代步引入虚拟项，冻结Poisson方程涉及的非线性项，然后利用仅含边界积分的基本解Trefftz 有限元模型进行求解.为建立基本解单元列式，需在单元域外设置若干虚拟的源点[31]，利用源点和单元节点的坐标信息计算基本解，进而构建域内插值函数.与完备解方法类似，虚拟源点的最优数目也需要通过数值算例验证.但Wang和Qin建议，采用与单元节点相同数目的虚源点能够保证求解精度.由于大多数工程和物理问题的基本解已获知，因此基于基本解构建域内插值函数的方法十分便捷.表2列出了基于完备解和基本解构造Trefftz插值函数的相关情况对比.另外，关于基本解Trefftz有限元法的综述文献可参见文献[32].
表2完备解和基本解Trefftz单元情况对比
对比项完备解Trefftz单元基本解Trefftz单元
齐次解获取较难.获取容易，因为大多数工程和物理问题的基本解已
获知.
Trefftz函数通过直接截断完备解系构建，但要保证阶次连续且满足截断准则.同时，需借助数值算例验证合理的
Trefftz项数.需按某种规则在单元域外围布置虚源点，以便计算基本解.虚源点个数一般取成与单元节点数相同.
零能模式完备解系含有与零能模式相关的刚体运动项，构造单元时需人工去除.计算单元域内场变量时再用最
小二乘法等方法恢复.基本解本身不含与零能模式相关的刚体运动项，故此构造单元时无需人工去除.计算单元域内场变量时，可用与完备解Trefftz单元相同的恢复方法.
2源项的无源化处理
基于前述两个独立的场变量插值模式，结合杂交泛函可推得Trefftz单元刚度方程.对于无源项的工程和物理问题，其控制方程是齐次的，在单元刚度方程推导过程中，采用高斯散度定理可直接消除杂交泛函涉及的域积分[2,12].然而，若源项(如弹性问题中的体力或位势问题中的源汇等)存在，问题的控制方程是非齐次的，高斯散度定理虽仍可消除泛函中的原有域积分，但源项还会衍生另一个新的域积分并出现在单
元节点载荷列阵中[2,33-35].以各向同性平面位势问题为例，其控制方程为Poisson方程，这里考虑第一类和第二类边界条件，则与该问题等价的杂交泛函Πe为
Πe=
1
2
Ωe
q21+q22
dΩ−
Γeu
q e¯u e dΓ+
Γeq
(¯q e−q e)˜u e dΓ−
ΓeI
q e˜u e dΓ(3)式中，q1和q2为沿x,y坐标方向的势流，Γeu和
436力学与实践2017年第39卷
Γeq分别为给定位势¯u e和法向势流¯q e的边界，ΓeI
为单元交边界，且有Γe=Γeu∪Γeq∪ΓeI，变量上方
的横线表示给定值.应用高斯散度定理，式(3)可改
写为
Πe=1
2
Γe
q e u e dΓ−
1
2
Ωe
bu e dΩ−
Γe q e˜u e dΓ+
Γeq
¯q e˜u e dΓ(4)
对比式(4)与式(3)可以看出，上述操作消去了原泛函右边第一项的域积分(式(3)右侧第1项)，但出现了与源项b相关的新域积分(式(4)右侧第2项).此时，单元域内插值函数N ej不再满足原问题的控制方程，须增加源项b诱发的特解部分 u e.众所周知，在传统有限元法中，由于单元列式中含有域积分，因此单元形状受雅可比矩阵控制，相邻单元边内角不能接近或等于180◦，即单元不能过分扭曲，否则会导致计算结果失真甚至无法获得解答.为避免出现传统有限元法面临的尴尬境地，也为保持Trefftz有限元法仅含边界积分的独特优势，研究者们提出了一些处理非齐次方程的方法，但共同作法是将非齐次问题转化为齐次
问题求解.目前，主要有格林函数法和径向基函数法.格林函数法解决了带常体力的平面弹性问题[2].为了处理任意形式的非齐次项情形，Qin[2]和Wang等[33]基于双重互易边界元法的概念提出了径向基函数法，将径向基函数定义为欧几里得距离的变量，来近似表示源项的函数分布.然后对控制方程求解析积分直接获得相关特解，这样原问题的求解过程就归为寻求特解和齐次解的问题.基于这种求解策略，Weiber 等[8,19,28,36-37]分析了一系列Poisson类方程的问题.王克用等[38-39]通过坐标变换和径向基函数分析了有源项正交各向异性平面和轴对称位势问题. Moldovan[40]构造了Trefftz应力元并分析了非齐次双曲边值问题.尽管径向基函数在处理非齐次问题时非常便捷，但它可能会造成病态的系数矩阵[41]，不过采取适当的规则化措施[42]避免解答严重失真.
3特殊功能单元
Trefftz有限元法最大的亮点之一是特殊功能单元.借助这种单元无需局部网格细分，只需在域内插值函数上作调整，就能捕捉各种奇异性或局部效应(如不连续载荷、角点、裂纹、夹杂和孔洞等)，凸显出传统有限元法望尘莫及的效率和优势[2,31-32].如前所述，杂交Trefftz单元的域内插值函数精确满足问题的控制方程.若域内插值函数还满足奇异或局部效应区域的边界条件，则可构造相应的特殊功能单元.在现有文献中，研究载荷奇异性的特殊功能单元的文章较少，较早的工作是由Jirousek等[43]完成的，他们推得平面应力状态下的特解构造了集中载荷单元和分布载荷单元.近期，Wang等[44]针对不连续载荷引起的局部效应构造了三种特殊功能单元，分别用于分析弹性体上的点载荷、线载荷和面载荷(图2)，均是基于
适当的局部基本解建立的.
在工程实际中，许多结构含有各种各样的孔洞，如螺孔、工艺孔等.较早研究孔洞单元的是Piltner[45]，他利用复变Laurent级数推导了椭圆孔单元，该单元能退化为圆孔和内裂纹单元，属于Tr-efftz型特殊功能单元的范畴，但当时Piltner并没有提出Trefftz有限元的概念.在Piltner[45]工作基础上，王克用[34]引入旋转映射函数构造了可任意调整椭圆孔倾角的特殊功能单元，并应用于接触问题分析中.Leconte等[46]探讨了8节点圆孔单元由线弹性问题推广至非线性问题(冲击)的可能性.Wang 等[47]分析了平面弹性问题，以格林函数作为域内插值函数构造了圆孔单元
.
(a)点载荷单元(b)线载荷单元(c)面载荷单元
图2不连续载荷特殊功能单元
第5
期王克用等：Trefftz 有限元法的研究进展437
若孔洞内填塞有与其形状相同、力学性质不同的材料就形成了夹杂，夹杂与基体之间存在相互作用，如颗粒/纤维增强复合材料等非均质材料.对于这种情形的Trefftz 有限元建模与分析，始于Zhang 等[48-50]的工作.他们先后提出了含图形或椭圆形弹性夹杂/刚性夹杂单元，其构造原理是将原单元分解为两个边值子问题，进而建立单元节点力与位移的关系(图3).观察发现，Zhang 等[48-50]的工作与Voronoi 单胞有限元法十分相似，而这种有限元法已经成为模拟非均质材料微观力学特性的强有力工具.张洪武等
[51]
基于参数变分原理构造了含夹
杂Voronoi 单元，分析了非均质材料中夹杂对其宏观等效弹塑性力学性能的影响.近期，Dong 等
[52]
在
Voronoi 单胞有限元理论
[53]
基础上，提出了Trefftz
型Voronoi 单胞有限元法，通过引入单元的特征长度，可以有效处理病态方程情形.在单元域内和边
界上，借助特征长度这个参数可将Trefftz 完备函数限定在[0,1]区间内，这样能确保系统方程组处于良态，无需额外的规则化技术进行处理就能轻松求解.实际上，单元特征长度所起的作用与杂交Trefftz 有限元法中通常采用的局部坐标系(o,x,y )[2,12,38-39,43]类同(图1).之后，Dong 等
[14]
将其工作推广至
含椭圆孔洞以及弹性/刚性夹杂单元.紧接着，Dong 等[15-16]还构造了三维Trefftz 多边形单元研究球形和椭球形孔洞、夹杂的非均质材料.Wang 等
[54]
则构造了特殊多边形Voronoi 纤维/基体单元，用于分析天然纤维复合材料的热效应.从构造原理来看，孔洞单元可视为同形状夹杂单元的特殊情形.此外，Cao 等
[55]
利用适当的局部基本解研究了含缺陷的平面压电问题，Wang 等[56-57]研究了平面弹性问题的特殊功能单元.关于基于基本解的Trefftz 特殊功能单元的综述可参见文献[58].图3非均质结构与Voronoi 单元
此外，Bishay 等[59]构造了一系列用于分析带有缺陷、孔洞以及弹性电介质/压电夹杂的压电材料单元.
每种单元的外边界条件应用变分原理、配点法或最小二乘法强制满足，而孔洞/夹杂外围应力/电荷自由边界条件则采用配点法/最小二乘法或者特殊解系得到满足.Hennuyer 等
[60]
构造了一种杂交Trefftz 超单元，用于模拟航天器结构在碰撞和冲击载荷作用下铆钉装配体的应力集中情况.为了考虑孔边带有初裂纹情形，他们增加了超单元的节点
数，尽管改进后的高阶超单元精度有所降低，但仍能对力场进行良好的描绘.Kunter 等
[61]
基于既
满足控制方程又满足裂纹边界条件的精确解，构造了Trefftz 裂纹单元，且推导了一个含有Dugdale 条状屈服区域的二维直裂纹的特解.Chen 等
[62]
基于
Hellinger--Reissner 变分原理提出了一种多边形角点
单元，它虽然可分析含有相互作用的双菱形、双正方形和双矩形孔洞的平面弹性问题，但其只是利用两个角点单元拼接成孔洞而已.这种单元能够较好地

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。