与平行计算的MATLAB中加速有限元分析

有限元方法是一种强大的数值技术，用于在一系列复杂的科学和工程应用中解决普通和部分微分方程，例如多域分析和结构工程。它涉及在构造之前将分析域分解成分立网格，然后求解在网格元件上构建的等式系统。随着网格精制的等式的数量增加，使得有限元方法计算非常密集。然而，该过程的各个阶段可以容易地并行化。

在本文中，我们执行耦合的耐静电微机电（MEMS）装置的电动机械有限元分析。我们将并行计算技术应用于机械分析阶段的最具计算密集的部分。使用40工人¹设置，我们将减少机械分析所需的时间，从近60小时到近60小时到少于2小时的大约一百万个DOF网。

其中一个最常见的MEMS开关是悬臂系列（图3）。这包括悬挂在接地电极上的光束。

图4显示了建模的几何形状。顶部电极的长度为150μm，厚度为2μm。杨氏模量E是170 GPA，泊松比υ为0.34。底部电极的长度为50μm，厚度为2μm，并且从顶部电极的最左端设置100μm。顶部和底部电极之间的间隙为2μm。

当在顶部电极和接地平面之间施加电压时，在导体的表面上诱导静电电荷，这导致作用在导体表面的静电力。由于接地平面是固定的，因此静电力仅变形顶部电极。当光束变形时，电荷在导体表面上重新分配。所得到的静电力和光束的变形也改变。该过程持续到达到均衡状态。

应用FEM以耦合电力学分析

为简单起见，我们将使用基于松弛的算法而不是牛顿方法来耦合静电和机械域²。步骤如下：

1.在可动电极上用恒定电位V0求解非成形几何形状中的静电FEA问题。

2.使用沿可动电极的电荷密度的计算值计算机械解决方案的负载和边界条件。

可动电极上的静电压力由\ [p = \ frac {1} {2 \ epsilon} | d | ^ 2 \]给出

在哪里，
\（| d | \）=电量密度的幅度
\（\ epsilon \）=可动电极旁边的电介电性

3.解决机械FEA来计算可动电极的变形。

4.使用可动电极的计算出位移，更新沿可动电极的电荷密度。

\ [\ left | d _ {\ mathrm {def}}（x）\ \ loct |\左\左| d_0（x）\ lex | \ frac {g} {g - v（x）} \]

在哪里，
\（| d _ {\ mathrm {def}}（x）| \）=变形电极中的电量密度的大小
\（| d_0（x）| \）=未变形电极中的电量密度的幅度
\（g \）=在没有致动的情况下可移动和固定电极之间的距离
\（v（x）\）=可动电极在其轴轴上的位移

5.重复步骤2 - 4，直到最后两个迭代中的电极变形值收敛。

静电分析（步骤1）涉及五个步骤：

• 绘制悬臂开关几何形状。
• 将Dirichlet边界条件指定为20V至可动电极的恒定电位。
• 网状外部域。
• 解决外部领域的未知电位。
• 绘制解决方案（图5）。

我们可以通过部分微分方程工具箱™接口快速执行这些任务中的每一个。

机械分析

机械分析涉及五个步骤：

• 绘制域名
• 衍生元素级方程
• 集会
• 施加边界条件
• 解决方程式

绘制域名

在此步骤中，我们将系统域分散到更小的域或元素。离散化允许我们代表域的几何形状，并近似于每个元素的解决方案，以更好地代表整个域上的解决方案。元素数量确定问题的大小。

衍生元素级方程

在该步骤中，我们假设在所选点或节点处的元素上的差分方程的近似解。在我们的示例中，解决方案是以\（x \）和\（y \）方向（2d分析）中的\（\ phi \）位移的离散值来确定的。节点处的未知主字段的数量称为该节点处的自由度（DOF）。

现在将控制微分方程应用于单个元素的域（图6）。在元素级别，控制方程式的解决方案被近似于元素域\（d ^ e \）上的连续函数替换，以未知的节点值表示（\phi_1 \），\（\ phi_2 \），\（\ phi_3 \）的解决方案\（\ phi \）。然后，可以为元素配制\（\ phi_1 \），\（\ phi_2 \）和\（\ phi_3 \）方程系统。

集会

通过组合每个元素的解决方案来获得系统的解决方程，以确保每个节点处的连续性。在我们的示例中，元素级刚度和力矩阵（\（k_e \）和\（f_e \））被组装以创建全局刚度矩阵（\（k））和力矩阵（\（f \））在整个域上

施加边界条件

我们强加了边界节点的必要边界条件，解决了全局方程式系统。解决问题的解决方案\（\ phi（x，y）\）成为分段近似，以\（\ phi \）的节点值表示。来自组装过程的线性代数方程的系统：\（k \ phi = f \）。

实际工程问题与包含数千方程的系统并不罕见。并行计算技术可以大大减少组装矩阵所需的时间，并计算出这种大规模的问题的解决方案。

解决方程式

我们解决了\（x \）和\（y \）方向上的位移的全局方程式\（k \ phi = f \）系统。

并行化问题

我们在MATLAB中实现了在有限元方法中定义的每个步骤^®功能。网格化完成使用mesh2d.，基于Matlab基于Mesh生成应用程序的2D几何形状。MATLAB中的分析工具表明，最耗时的操作属于刚度矩阵组装步骤。

此步骤涉及每个元素的三个主要操作：
1.从元素级解决方程计算元素刚度矩阵（\（k_e \））。\（k_e \）是大小的\（n_e \ times n_e \）
在哪里，

2.将\（k_e \）矩阵值的本地位置映射到全局刚度矩阵中的位置。
3.使用具有元素刚度矩阵值的地图填充全局刚度矩阵（\（k））。

表1.不同DOF的组装时间和总执行时间的比较。

随着元素的数量增加，迭代的数量和\（k \）的大小。图7和表1将组装到具有不同DOF系统的串行模式的总执行时间所花费的时间。显然，大会是最耗时的部分，占据系统具有高DOF的总执行时间的99％以上。

刚度矩阵组件

幸运的是，最终组装的\（k）矩阵与在循环中选择元素的顺序无关。我们可以通过在多个MATLAB工人跨多个MATLAB工人分发计算来评估几个元素与全局刚度矩阵（\（k））的贡献。装配操作通常在串行中执行为了-Loop，通过每个元素步骤并确定其对全局刚度矩阵的贡献。我们只是转换序列为了- 并行为了- 使用它议案并行计算工具箱™构造（图8）。

在我们的示例中，全局刚度矩阵是元素刚度矩阵对循环的所有迭代的贡献之和。这议案构造使我们能够自动处理这种缩减分配（通常是表单\（r = f（\ f（\ text {expr}，r））））。

为了展示通过使用并行方法实现的性能增益，将问题从粗地网格缩放到超级格式网格。粗滤网包含大约128个元素，共150分。我们改进了网状物，直到它包含856,800个元素，861,122 DOF。当我们改进网格时，悬臂梁的自由端的位移会聚

对于并行方法，我们使用带有一个头部节点和5台机器的计算机群集，每个计算机都有以下配置：双英特尔^®Xeon.^®1.6 GHz四核处理器（每台机器8芯，共有40个核心），13 GB RAM带窗户^®64位操作系统。每台机器都跑了8名MATLAB工人，共有40名工人。要测量串行执行时间，我们使用了在头节点上运行的单个MATLAB工作台。使用64位OS使我们能够创建大稀疏矩阵（最多861,122 x 861,122），而不耗尽内存限制。在大多数有限元应用中，所得到的k在自然界中稀疏。

图9A和9B在串行（RED）和并行（绿色）执行模式之间的DOF中的增加，将总执行时间进行比较。表2总结了结果。

表2.图9A和9B中使用的代表性数据。

请注意，对于具有少量DOF的系统，与这些操作的执行时间相比，分发操作的成本要高得多。因此，当系统仅有250个DOF时，串行执行实际上比并行执行更快。

图9B显示了最多约400个DOF，两个曲线相交的点，串行模式执行比并行模式更快。我们通过在此之后仅切换到并行模式，我们会看到性能改进。实际执行时间和交叉点取决于若干因素，包括所涉及的MATLAB函数的执行速度，工作人员的处理速度，网络速度和工人数量，可用内存和系统负载。

概括

在本文中，我们展示了一种简单的方法来并行化FEA应用程序。我们开始通过分析串行代码性能，专注于我们设置的最具计算密集的部分。通过简单的代码，我们能够显着提高应用程序性能，减少时间从60小时从60小时分析到40小时到少于2小时的时间。

¹工人是Matlab计算引擎，与Matlab会话分开运行。

²P.S.SUMANT，N.R.Aluru和A.c. Cangellaris，“静电动作MEMS的快速有限元建模的方法。”工程中的数值杂志2009;77：1789-1808。