英语原文共 12 页，剩余内容已隐藏，支付完成后下载完整资料
应用数学与计算
基于边界单元的2D和3D有限元矩阵的MARLAB快速组装
文章信息：
关键词：MATALAB代码向量化、有限元方法、边界单元法、Raviart-Thoma单元。
摘要：
对于如何组装MATLAB中的有限元刚度矩阵和质量矩阵，我们提出了一个有灵活、有效的方法。我们将这种方法应用在边界有限元离散化的问题中。典型的边界有限元是用于H（div）空间中离散化的Raviart-Thomas单元和H（curl）空间中离散化的Nedelec单元。我们解释了向量化思维，并且对一种可用程度十分高的MATLAB代码进行了测评，这种代码十分便捷，并且可以拓展。
1、介绍
在H Sobolev空间中，椭圆问题以弱形式的方式表现出来的，并且用节点(nodal)有限元函数将其离散化。椭圆问题包含了标量和向量参数的全部梯度算子。Rahman和Valdman在论文[11]中解释了线性节点有限元刚度矩阵在MATLAB中组装路线的有效向量化。这篇论文的重点是将这种方法推广到包括高阶单元在内的任意的有限元和需要用散度算子和旋转度算子处理向量问题。这种问题出现在电磁领域，也和机械领域中各种多重或者二重问题有关。这些问题的弱形式定义在H(div)和H(curl) Sobolev空间中。有限元的离散化是根据边界单元完成的，典型的这些单元有H(div)问题中的Raviart-Thomas单元和H(curl)问题中的Nedelec单元。边界单元基函数并不是定义在二维三角形或三维四面体的网格上，而是在边和面上。边界单元仅仅在单元的边界上才具有连续性——H(div)问题的法向向量分量连续性和H(curl)问题的切向向量分量连续性。
目前，用于处理H(div)和H(curl)问题的有限元方法和它的工具库都很好地记录了下来，比如例子[16]，通过分层依据对高阶多项式进行了定义。用户会发现许多软件代码其实都是用面向对象语言编辑的（比如NGSOLVE或者HERMES），这种语言允许用曲面边界单元定义高阶单元。这种语言在处理高度复杂运算方面是十分出色的，并且通过用户界面提供确定的自由度。然而，理解或修改这种代码是十分困难的，除非用户十分熟悉这种代码。我们相信，MATLAB代码对那些希望熟悉边界单元法，并且有它们自己的实现方法的学生和调查员是十分方便的。我们认为，最低级的线性边界单元仅仅定义在二维三角形和三维四面体上。然而，将这种代码延伸到高阶单元也是很简单的，因为不论单元是高阶还是低阶，组装路线是几乎一样的。
有大量的论文旨在实现MATLAB节点单元组装路线的向量化。在论文[8]中，作者也讨论了Raciart-Thomas单元在三维空间的状况，但是没有提供程序代码。iFEM库中有用于实现不同的线性和高阶单元组装路线的有效方法。论文[2]考虑了Raviart-Thomas单元法（以一种没有向量的方式），这促进了多重网格求解器在H(div)强函数最小化方面的应用（经由第二作者）。
我们的实现方法推广了论文[11]中的方法，使该方法能够处理任意的仿射有限元。它也是基于MATLAB长向量和阵列算子的，并且它十分适合处理大规模问题。在一台有着不错的处理器和足够系统内存的电脑上，有限元矩阵的二维或三维组装是很快的。例如，一个大约1千万行的二维矩阵组装花费了不到1分钟。矢量计算常常需要更多的系统内存，但是只有在系统内存变满的情况下，性能才会下降。这篇论文中提到的软件可以在http://www.mathwork.com/matlabcentral/fileexchange/46645免费下载。该网址还有一个二维和三维线性节点有限元的工具，以及它的推广方法，因此有兴趣的读者可以比较一下论文[11]中不同代码的性能。我们广义方法的关键思路是将仿射模型值函数的标量和向量的积分过程向量化。这也使快速评估范函数成为可能。
论文分为以下几个部分：在第二章中，我们简要介绍了如何实现线性边界单元法。第三章中，我们把有关单元和向量工具细节方面的特殊结构完成了。我们也展示了向量组装路线在时间和拓展性方面的表现。第四章介绍了两个边界有限元的应用：后验误差分析的强函数的最小化和电磁问题的解决方法。
2、线性边界单元
我们用Omega;表示一个开放的、有边界的、连通的Rd利普希茨域，其中，disin;{2,3}，表示了空间的维度。向量值函数w:Omega;-Rd的散度（二维和三维）和旋度（三维）被定义如下：
我们考虑了两类二维旋度算子——向量旋度算子和标量旋度算子。
适用于一个标量函数f：Omega;-R和一个向量函数w：Omega;-R2.在文献中，旋度算子常被成为“同步”，并且常被表示为▽。这些算子会产生Sobolev空间：
L2表示利普希茨积分函数的平方。我们将用||.||：=||.||L(Omega;)表示标量和向量函数的L2规范。假设Omega;被三角形（二维）和四面体（三维）网格划分，那么Raviart-Thomas和Nedelec单元代表H(div)和H(curl)空间的基函数。
在最低阶（线性）Raviart-Thomas和Nedelec单元情况下，有一个全局自由度，即，一个网格的每个边或者面都相关的全局基函数。基于这个解释，可以想到全局Raviart-Thomas基函数和二维Nedelec基函数仅有两个单元是非零的，这两个单元共享与基函数相关的边或面。在三维中，Nedelec全局基函数的所有单元都是非零的，这些单元共享相关的边，并且它们的数量通常超过两个。
我们用eta;RT、eta;Ned和Omega;空间变量X=(x1,x2,x3)T表示全局边（面）基函数。参考基函数的符号和空间变量可以通过简化位的求和获得，即，X表示参考单元K的空间变量。我们将用二维单元三角形和三维单元四面体作为参考单元。我们用ei表示参考三角形或四面体的边，用fi表示参考四面体的面。边和面的编号以及自由度的编号可以在图1中看到。接下来的Fk表示从参考单元K到网格单元K的仿射单元分布像：Fk(x):=Bkx bk。
图1：参考位形K上线性边界单元的自由度
图2：由节点和边组成的单元，即，二维线性有限元的全局自由度编号
通过一个三元集合{K,R,A}定义一个有限元，其中，K表示参考位形，R表示定义在参考位形的函数的有限空间，A表示线性无关自由度的集。我们已经选择了参考位形因此我们还需要一个映射，该映射接受来自R的函数，并将它们从K映射到网格上的单元K，这种映射称为Piola映射。
2.1、Raviart-Thomas单元
线性Raviart-Thomas单元基于空间(see.e.g,[9,12])和自由度uisin;R
每个二维的边为e，三维的面为f，和相应的参考单元为K。这种情况下，二维空间有3个自由度，三维空间有4个自由度。这里的ni表示边和面的法向单位向量。这里必须选择两个单位法向量中的一个来使用。外部单元法线的选择标准如图1所示。alpha;规则给予了我们Raviart-Thomas单元的参考基函数。
为了保护相关基函数的法向连续性，我们需要使用所谓的Piola映射。值和分散值的映射如下所示：
2.2Nedelec单元
线性Nedelec单元基于e.g[9,13]空间，
并且Uisin;R的自由度和边界单元有关。
对于相关位形K中的每个边e。二维中有三个自由度，三维中有六个。T是e的切向单元向量。类推到Raviart-Thomas单元中，需要确定所用切向单元向量的方向。我们的选项列举在图1中。alpha;规则给予了我们Nedelec单元的相关基函数。
同样，我们需要用一个Piola映射来确保切向连续性。值的映射如下：
旋度的映射取决于维度。
2.3局部自由度方向。
为了获取eta;的全局基函数，前面章节介绍的变换是不够的。一个全局基函数至少与两个变量有关。没有相关文献能确认这些不同单元的局部自由度方向是相同的。为了让Raviart-Thomas和Nedelec单元产生相同的函数，这些方向必须相同，而这些函数的切向分量和法向分量在单元表面是连续的。
以二维Raviart-Thomas单元为例，Km和 Kn是共享一条边的两个单元，eta;RT是这条边的全局参考基函数。为了得到全局基函数，我们用Km和 Kn表示由K变换到Km和 Kn的参考基函数。
通过观察自由度，我们发现我们总是在用外部单位法线来计算局部基函数。如果我们只是用转换(2)，边缘值的法向分量或许会跟变得彼此都不一样。这取决于单元分布像Fkn和Fkm是否保持方向不变。如果行列式BKN大于零，且行列式Bkm小于零，那么单元分布像Fkn就一直是参考基函数的逆时针方向，并且Fkm是顺时针方向。这意味着在公共边上，方向是相同的，转换(2)相对两个单元也是足够的。然而，如果情况相反，那么公共边的方向将是相反的，并且其中的一个转换必须乘以-1.因此全局基函数要通过如下方式获得：
图3：当两个单元都定为逆时针且cgt;a时，二维边界单元的方向沿着一条边。粗线表示正向
我们将这些值成为符号数据，它们与两个单元的中每一个都有关。注意，综上所述，全局基函数可以被简化为：
但是我们将使用公式(7)，因为它更容易用程序代码实现。
这个情况同样适用于三维Raviart-Thomas单元和二维Nedelec单元。至于三维Nedelec单元，用户需要更加小心，因为在一个规模相对较大的单元中，全局基函数可能是非零的，并且参考方向和边有关。还要注意，必须使用相同的符号数据来变换基函数的旋度或散度。
2.4有限元矩阵
我们对质量矩阵MRT，MNeD和刚度矩阵KRT和KNeD的组装十分感兴趣。它们定义如下：
指数i和j是自由度的全局编号，即，它们跟一个网格的边或者面有关。通过使用Piola映射（需要确定正确的方向），我们能够使用参考单元来组装局部矩阵。
通过使用公式(7)，与全局矩阵MRT和KRT有关的局部矩阵可以用单元K来进行计算(Kisin;T):
K是参考单元。指数k和l运行与所有的局部基函数，在二维中，k,lisin;{1,2,3}，在三维中，k,l属于{1,2,3,4}。
同样地，通过使用公式(4)，同时也要考虑正确的方向（参见2.3部分），与整体质量矩阵有关的单元质量矩阵MNed可以用单元K计算，K属于T：
通过使用公式(5)和(6)，和整体刚度矩阵有关的单元刚度矩阵KNed可以通过如下方法计算：
指数k和l适用于所有的局部基函数，在二维中，k,lisin;{1,2,3}，在三维中，k,l属于{1,2,3,4}。
3.边界单元的实现
我们用#w表示集合w的单元数，相对的，用N，ε，rho;和T表示点，边，面和单元。注意，面rho;只存在于三维。为了实现边界单元，我们需要用下列结构表示网格。第二列说明了结构的大小，第三列说明了结构的意义。
点2坐标 #Ntimes;2/3 二/三维中的两/三个坐标决定点
边2节点 #εtimes;2 二/三维中的两个节点决定一条边
面2节点 #Ftimes;3 三维中的3个点决定一个面
通过这些可用的矩阵，我们能用列表中的点，线或面来描述每一个单元。
单元2节点 #Ttimes;3/4 二/三维中的3/4个节点决定单元
单元2边 #Ttimes;3/6 二/三维中的3/6条边决定单元
单元2面 #Ttimes;4 三维中的四个面决定单元
在二维中，线性Raviart-Thomas单元和Nedelec单元都有一个自由度，这个自由度和每个参考三角形的三条边有关，一共三个自由度。在三维中，线性Nedelec单元有一个和四面体的六条边相关的自由度，并且Raviart-Thomas单元将有一个和四面体的四个面相关的自由度。因此，全局自由度的编号由行指数 [边2的点]或[面3的点]确定。然后，对于网格的特殊单元K，跟它有关的全局自由度由[单元2边]或[单元2面]确定。全局自由度编号的节点单
全文共11325字，剩余内容已隐藏，支付完成后下载完整资料
资料编号：[11622]，资料为PDF文档或Word文档，PDF文档可免费转换为Word
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
