用并行计算加速MATLAB有限元分析

有限元法是一种强大的数值技术，用于求解常微分方程和偏微分方程在一系列复杂的科学和工程应用中，如多域分析和结构工程。它涉及到在构建和求解基于网格元素的方程组之前，将分析域分解为离散网格。随着网格的细化，所涉及的方程的数量也在增加，这使得有限元方法的计算量非常大。但是，该过程的各个阶段可以很容易地并行化。

在这篇文章中，我们对一个静电驱动的微机电机械(MEMS)装置进行了耦合的机电有限元分析。我们将并行计算技术应用于力学分析阶段中计算量最大的部分。雇用40名员工¹安装完成后，我们将把大约100万自由度网格的力学分析所需的时间从近60小时缩短到不到2小时。

最常见的MEMS开关之一是悬臂系列(图3)。它由悬挂在接地电极上的梁组成。

图4显示了建模的几何图形。顶部电极长150μm，厚2μm。杨氏模量E为170 GPa，泊松比υ为0.34。底部电极长50μm，厚2μm，距顶部电极最左端100μm。上下电极间距为2μm。

当在顶部电极和接地面之间施加电压时，导体表面会产生静电荷，从而产生与导体表面正常作用的静电力。由于地平面是固定的，静电力只使顶部电极变形。当光束变形时，电荷在导体表面重新分布。所产生的静电力和梁的变形也会发生变化。这个过程一直持续到达到平衡状态。

有限元在机电耦合分析中的应用

为了简单起见，我们将使用基于松弛的算法而不是牛顿方法来耦合静电和力学域²．具体步骤如下:

1.求解可动电极上电位V0不变几何形状下的静电有限元分析问题。

2.利用沿活动电极电荷密度的计算值，计算力学解的负载和边界条件。

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

在那里,
\(|D|\) =电通量密度的大小
\(\epsilon\) =活动电极旁边的电介电常数

3.通过力学有限元分析计算可动电极的变形。

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

\[\左| D_ {\ mathrm {def}} (x) \右| \大约\左|数(x) \右| \压裂{G} {G - v (x)} \]

在那里,
\(|D_{\ mathm {def}}(x)|\) =变形电极中的电通量密度的大小
\(|D_0(x)|\) =未变形电极中的电通量密度的大小
\(G\) =无驱动时，活动电极与固定电极之间的距离
\(v(x)\) =活动电极在x位置沿其轴的位移

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

静电分析(第一步)包括五个步骤:

• 绘制悬臂开关几何形状。
• 指定狄利克雷边界条件为对可移动电极的恒定电位为20v。
• 网格外部域。
• 求解外部域的未知电势。
• 绘制解决方案(图5)。

我们可以通过偏微分方程工具箱™界面快速执行这些任务。

粒度分析

力学分析包括五个步骤:

• 域网格化
• 推导元素级方程
• 组装
• 施加边界条件
• 解方程

领域网格划分

在这一步中，我们将系统域离散为更小的域或元素。离散化让我们表示域的几何形状，并近似每个元素的解，以更好地表示整个域的解。元素的数量决定了问题的大小。

推导元素级方程

在这一步中，我们假设在选定的点或节点上的微分方程有一个近似解。在我们的例子中，解是根据\(\phi\)位移在\(x\)和\(y\)方向上的离散值确定的(2D分析)。节点上未知主字段的数量称为该节点的自由度(DOF)。

控制微分方程现在应用于单个元素的域(图6)。在元素级别，控制方程的解被一个连续函数取代，它近似于在元素域(D^e)上的\(\phi\)分布，用解\(\phi\)的未知节点值\(\phi_1\)、\(\phi_2\)和\(\phi_3\)表示。这样就可以用\(\phi_1\)、\(\phi_2\)和\(\phi_3\)来表示该元素的方程组。

组装

我们通过组合各元素的解方程来保证每个节点的连续性，从而得到系统的解方程。在我们的例子中，元素级刚度和力矩阵(\(K_e\)和\(F_e\)被组装起来，在整个域上创建一个全局刚度矩阵(\(K\))和一个力矩阵(\(F\))

施加边界条件

在边界节点处施加必要的边界条件，求解全局方程组。该问题的解\(\phi(x,y)\)变成了分段逼近，用\(\phi\)的节点值表示。一个线性代数方程组的结果由装配程序:\(K \phi = F\)。

在实际工程问题中，包含数千个方程的系统是很常见的。并行计算技术可以大大减少组装矩阵和计算如此大规模问题的解决方案所需的时间。

解方程

我们求解了位移在x方向和y方向上的全局方程组K = F。

并行解决问题

我们在MATLAB中实现有限元法中定义的每个步骤^®函数。网格划分是用mesh2d，一个基于MATLAB的二维几何网格生成应用程序。MATLAB中的Profiling工具表明，最耗时的操作属于刚度矩阵装配步骤。

这一步涉及到对每个元素的三个主要操作:
1.根据单元级解方程计算单元刚度矩阵(\(K_e\))。\(K_e\)的大小\(N_e \乘以N_e\)
在那里,

2.将\(K_e\)矩阵值的局部位置映射到它们在全局刚度矩阵中的位置。
3.使用带有元素刚度矩阵值的映射填充全局刚度矩阵(\(K\))。

表1。不同自由度装配时间和总执行时间的比较。

随着元素数量的增加，迭代次数和\(K\)的大小也会增加。图7和表1比较了不同DOF系统在串行模式下的装配时间与总执行时间。显然，装配是最耗时的部分，当系统具有高DOF时，装配占总执行时间的99%以上。

刚度矩阵组合

幸运的是，最终组装的\(K\)矩阵与循环中元素选择的顺序无关。我们可以通过将计算分布到多个MATLAB工作者，同时评估几个元素对全局刚度矩阵(\(K\))的贡献。装配操作通常以串行方式执行为-loop，它遍历每个元素并确定其对全局刚度矩阵的贡献。我们只需转换串行为-loop到并行为-循环使用parfor在并行计算工具箱™中构造(图8)。

在我们的例子中，全局刚度矩阵是在循环的所有迭代中元素刚度矩阵的贡献之和。的parfor构造使我们能够自动处理这种约简赋值(通常是\(r = f(\text{expr}， r)\)的形式)。

为了演示使用并行方法所获得的性能增益，将问题从粗网格扩展到超细化网格。粗网格包含约128个元素，总自由度为150。我们将网格细化到包含856,800个元素，DOF为861,122。当我们细化网格时，悬臂梁自由端位移收敛

对于并行方法，我们使用具有一个头节点和5台机器的计算机集群，每台机器具有以下配置:Dual Intel^®至强^®1.6 GHz四核处理器(每台机器8核共40核)，13gb内存^®64位操作系统。每台机器运行8个MATLAB工人，共40个工人。为了测量串行执行时间，我们使用了一个运行在头节点上的MATLAB工作器。使用64位操作系统使我们能够创建大型稀疏矩阵(高达861,122 x 861,122)，而不会遇到内存限制。在大多数有限元应用中，得到的K本质上是稀疏的。

图9a和9b比较了随着串行(红色)和并行(绿色)执行模式之间DOF的增加，总执行时间(以秒为单位)。表2总结了结果。

表2。图9a和9b中使用的代表性数据。

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

图9b显示了大约400 DOF，即两条曲线相交的点，串行模式执行速度比并行模式快。只有在这之后切换到并行模式，我们才能看到性能的提高。实际执行时间和交叉点取决于几个因素，包括所涉及的MATLAB函数的执行速度、worker机器的处理速度、网络速度和worker数量、可用内存和系统负载。

总结

在本文中，我们演示了并行化FEA应用程序的简单方法。我们首先分析了串行代码的性能，重点放在我们的设置中计算量最大的部分。通过简单的代码更改，我们能够显著提高应用程序的性能，将分析一个800,000DOF系统的时间从60小时缩短到40人设置的不到2小时。

¹worker是独立于MATLAB会话运行的MATLAB计算引擎。

²P.S. Sumant, N.R. Aluru和A.C. Cangellaris，“静电驱动MEMS的快速有限元建模方法”。国际工程数值方法杂志2009;77:1789 - 1808。