英语原文共 22 页，剩余内容已隐藏，支付完成后下载完整资料
用于SIPG线性弹性的快速原生 - MATLAB刚度组件
摘要
当在MATLAB中编写时，可以快速实施有限元法（FEM）并且与编译代码相比，行数更少。MATLAB也是一个有吸引力的环境，因为在它生成定制的科学计算例程中含有易于访问的内置函数库、有效的调试工具和一种生成脚本的简单语法。然而，人们普遍认为MATLAB对大问题的分析效率太低。在此基础上，提出了一种基于矢量化和闭塞算法来计算对称内罚不连续的Galerkin (SIPG) FEM的总体刚度矩阵算法。SIPG的总体刚度矩阵和传统的连续Galerkin近似计算的主要区别是要求对元素间的面项进行评估，这大大增加了计算工作量。本文着重研究了在计算时间上的面积分，并没有在现有文献中进行讨论。与文献中现有的优化的有限元素算法不同，本文只使用了原生的MATLAB功能，并且与GNU Octave兼容。该算法主要描述了基于齐次单元类型和多项式阶的网格的二维分析。同样的结构也适用于3D分析。对于106自由度大小的问题，局部刚度矩阵的2D计算速度至少快asymp;24倍，矢量化提高13.7倍，阻塞提高1.8倍。阻塞和矢量化的速度取决于计算机体系结构，本文将介绍两种体系结构的潜在改进范围。
1. 介绍
有限元分析（FEA）通常被工程师、数学家和科学家用作解偏微分方程的技术。MATLAB环境及其功能和调试程序库允许用少量行快速生成定制的FEA例程。例子包括Coombs等人[1]、西格蒙德[2]等人。然而，未优化的MATLAB脚本通常比未优化的编译代码慢得多[3]。本文演示了如何使用原生MATLAB生成FEA程序的优势，如果以对称的内部惩罚方法（SIPG）的优化形式编写，运行时间不会受到惩罚。
Dabrowksi等人[3]于2008年在MATLAB中优化FEA程序方面取得了重大进展。作者提出了MILAMIN，一种开源优化非本地MATLAB实现的连续Galerkin（CG）有限元分析代码，能够在一分钟内建立、解决和后处理具有106个自由度的2D非结构化网格问题。计算局部单元刚度的一种常用方法是通过一系列小矩阵乘法轮流计算每个局部单元矩阵。在MATLAB中创建MILAMIN算法时，作者认识到这种方法存在两个重大瓶颈。
首先，需要两个嵌套for循环来生成网格中的所有单元刚度矩阵。外循环，循环通过所有元素和内部循环，遍历所有高斯点。由于MATLAB循环本质上很慢，在计算大网格（超过106自由度）的刚度矩阵时，单元循环的迭代次数很大，这被认为是第一个瓶颈。其次，MATLAB中的矩阵计算由线性代数执行包（LAPACK），它调用基本线性代数子程序（BLAS）包。每次调用BLAS函数它的开销远远大于计算时间。在传统的FEA代码中，这个开销很大由于发生了许多小矩阵计算。
Dabrowsk等人[3]通过设计一个算法，这个算法可以同时为所有元素计算局部刚度矩阵中的条目，从而消除了两个瓶颈。这使得MATLAB for循环的迭代大小与问题的大小无关，并且也大大减少了BLAS调用的数量，从而减少了它们各自的开销。通过使缓存重用最大化，MILAMIN例程得到了进一步改进，这种技术被称为阻塞。这项工作已经通过在文献[4]中引入平行矢量化刚度矩阵计算而得到了扩展。
最近Rahman和Valdman[5]为线性节点形函数的元素的体积积分编写了一个快速的MATLAB脚本。重点是从具有标准有限元结构的非矢量化代码开始，然后通过矢量化提高其计算速度。其中一个关键特性是保留了代码的原始结构，这确保了可读性不会丢失，而代码优化通常就是这种情况。Dabrowksi等人[3]也强调了优化代码中缺乏可读性。此外，Anjam和Valdman[6]为典型的Raviart-Thomas元素制作了一个矢量化的MATLAB脚本，用于在H(curl)空间离散化的H(div)空间和Nedelec元素的离散化。Andreassen等人[7]提供了不同向量计算语言之间计算性能的比较和讨论，以组装FE总体刚度矩阵。Cuvelier等人[8,9]提出了一种更一般的多矢量语言矢量化例程。
在一个有限元分析代码中，一旦所有的局部单元刚度矩阵被计算出来，它们就组合在一起，作为单元拓扑的函数，形成一个稀疏的总体刚度矩阵。在本地MATLAB中，这是通过使用命令稀疏来实现的，该命令从三元组数据生成稀疏矩阵：行位置，列位置和相关联的值。本地表现很慢，Dabrowski等人 [3]使用sparse2一个非本地的MATLAB命令。其他稀疏矩阵MATLAB命令也已经在[10]中创建、研究、改进和讨论过，他们还提供了自己的改进和图形处理单元（GPU）实施。
本文实现了线性弹性的SIPG方法。不连续的Galerkin（DG）方法首先由Reed等人[11]引入，用于求解中子输运方程。Richter [12]提出将原始DG方法扩展到椭圆问题，包括线性对流扩散项。然而，不连续近似仅适用于对流项，对于二阶椭圆算子使用混合方法。Bassi和Rebay [13]介绍了对流和二阶椭圆算子的完全不连续逼近。
DG方法的一个突出特点是自由度是特定于元素的，允许在元素接口处进行简单的通信。特别是，hp-refinement由于其能够在元素接口处结合挂起节点而得到简化。这些特性使得DG方法非常适合高效的自适应求精以实现高保真模拟[14]。允许这种灵活性的代价是，与CG方法相比，对于相同数量和类型的元素，弱结合和自由度的积分项的数量更高。附加积分是面元连通性刚度项，它耦合元素之间的非共享自由度。这增加了产生总体刚度矩阵K[15]所需的计算次数，因此即使对于相对较小的问题，也需要有效生成K矩阵。
本文扩展了Dabrowski等人[3]提出的算法，包括SIPG的面部术语的优化整合，[16]，用于矢量化阻塞形式的线性弹性问题。这里所有的算法都是为本地MATLAB功能而设计的，这与文献[3-5]中大多数优化的MATLAB算法完全不同。Frank等人[17]在DG方法上唯一其他已知的矢量化，非阻塞的MATLAB代码。文献[17]中的作者将时间依赖扩散方程作为他们的模型问题，将其应用于2D局部DG公式中。这里我们在本地MATLAB中设计一个块矢量化代码，它利用SIPG中的对称性来模拟2D和3D中的线性弹性问题。
本文首先简要介绍了用于线性弹性的双线性SIPG公式。紧接着，将双线性SIPG格式重新组合为矩阵形式，可以在第2节中的向量化算法中计算。用于计算SIPG面刚度项的向量化算法在第3节中给出并讨论，体积积分被省略，因为它在[3]中已经被完全覆盖。用于计算SIPG面部刚度项的向量化算法对应于可从[18]中获得的Linear2D_DG.m脚本。计时结果，验证和讨论将在第4节中介绍，然后在第5节中作出结论。
2. 优化DG方法
2.1. 线性弹性问题的SIPG弱形式
在这里，我们考虑在Rd , d isin; {2, 3}上的有界Lipschitz多边形/多面体域Omega;上的下面的模型问题，边界part;ΩNcup;part;ΩD = part;Ω，其中part;ΩD和part;ΩN边界分别应用齐次Dirichlet和Neumann边界条件。对于小应变超弹性，问题的强大形式被定义为
gD和gN分别是应用的狄利克雷和诺伊曼边界条件。柯西应力张量定义为sigma; = part; (ε)/part;ε(u)，其中是超弹性的自由能函数，ε是小应变，u是位移，n是边界的法向单位向量。柯西应力张量也可描述为sigma; = Dε(u)，其中D是与应力和应变有关的材料刚度张量。
本文仅提供2D优化代码的描述，因此将省略3D元素空间和相应网格的描述。多边形有限元网格T在单元类型上是均匀的，并且通常是非结构化的。这里定义了两种元素类型，即三角形和四边形，但由于网格中只存在一种元素类型，因此这两种类型都称为K.多边形网格T由元素K组成，元素K可以是参考三角形的仿射图像元素映射FK：→K或双线性元素映射BK：→K下的四边形元素。三角形元素的均匀不连续Galerkin有限元空间定义为Wp(T ) = {w isin; [L 2 (Ω)] d: forall;K isin; T , w|K isin; Pp(K)},对于四边形单元Wp(T ) = {w isin; [L 2 (Ω)] d: forall;K isin; T , w|K isin; Qp(K)},，其中Pp(K)是小于或等于p的K上多项式的空间，Qp(K)是K上小于或等于p的多项式的空间。
我们用F(K)表示三角形的三个元素面的集合，或者表示元素K的四边形的四个元素面的集合。如果两个元素的交集F = part;cap; part;，,isin;T是一个段，我们称F为T的内部面。所有内部面的集合由FI（T）表示。类似地，如果元素Kisin;T和part;Ω的交点F = part;K cap; part;Ω是一个分段，我们称F为T的边界面。
用于近似模型问题（1）的SIPG方法现在以双线性形式引入，其中 part;ΩD上的齐次Dirichlet边界条件被强有力地应用并且不存在体力。找到位移uh isin; Wp(T )使得对于所有w isin; Wp(T )，有a(uh, w) = l(w)，其中
且
beta;是[19]中定义的线性弹性SIPG的惩罚项，hF是元素面的这个尺寸，并且
其中交点F isin; FI(T )上的和的元素面分别被称为F 和F-，其中n 和n-分别作为它们各自的向外法线。为方便起见(·, ·)和lang;·, ·rang;，其中（a，b）Omega;=int;Omega;ab和lang;a, brang;part;Ω = int; part;Ω ab.
2.2. SIPG方法的矩阵形式
既然已经描述了该问题的弱形式，可以将（2）中的应力，应变和位移表示为节点位移，形函数及其导数和材料刚度的函数。一旦表达出来，双线性形式中的每一项都可以重新表达为一组矩阵乘法，它们可以用来计算SIPG的刚度矩阵。实现矩阵形式的第一步是将元素位移uh分解成元素形函数Nn及其对应系数un的矩阵，使得uh = Nnun其中
un = [u1, v1, . . . , unen, vnen] T，nen是元素节点的个数，u和v分别是笛卡尔坐标系的x和y方向上的位移。类似地，测试函数可以表示为w = Nnwn。
由于小应变张量是uh的函数，应变也可以表示为一组矩阵乘法ε= LNnun，其中附加项
作为小应变偏微分算子[15]。从超弹性角度来看，柯西应力简单表示为sigma;=Dε= DLNnun，其中D是平面应力或应变刚度矩阵。（2）中的应力，应变和位移的矩阵形式，并设置
给出，
Neumann边界条件（3）中的测试函数项也表示为矩阵乘法
双线性形式中的每一项现在可以通过设定（10）等于（9），将括号相乘并除以wn得到
其中（KK <strong
全文共23796字，剩余内容已隐藏，支付完成后下载完整资料</strong
资料编号：[11184]，资料为PDF文档或Word文档，PDF文档可免费转换为Word
以上是毕业论文外文翻译，课题毕业论文、任务书、文献综述、开题报告、程序设计、图纸设计等资料可联系客服协助查找。
