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

有限元法是一种强大的数值技术,解决普通和偏微分方程在一系列复杂的科学和工程应用,如多域分析和结构工程。它包括将分析域分解为一个离散网格构建然后解决方程组建立在网格元素。方程所涉及的数量随着网格细化,使有限元法计算非常密集。然而,各个阶段的过程可以很容易地并行化。

在本文中,我们执行机电耦合有限元分析的静电驱动微机电(MEMS)设备。我们将并行计算技术应用到最计算密集型机械分析阶段的一部分。使用40-worker¹设置,我们将减少所花费的时间大约一百万个自由度的机械分析网格从近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 = \压裂{1}{2 \ε}| D | ^ 2 \]

在那里,
\ (D | | \) =电通量密度的大小
\ε(\ \)=电介电常数可动电极

3所示。解决机械有限元分析计算的变形可动电极。

4所示。使用可移动电极的位移计算,更新可动电极的电荷密度。

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

在那里,
\ (| D_ {\ mathrm {def}} (x) | \) =异形电极的电通量密度的大小
\(|数(x) | \) =未变形的电极的电通量密度的大小
\ \ (G) =移动和固定电极之间的距离没有驱动
\ (v (x) \) =位移可动电极的位置x轴

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

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

• 画出几何悬臂式开关。
• 狄利克雷边界条件指定为一个常数20 V移动电极的潜力。
• 网外部域。
• 解出未知电势在外部域。
• 情节的解决方案(图5)。

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

粒度分析

机械分析包括五个步骤:

• 啮合域
• 推导级方程式
• 组装
• 施加边界条件
• 解的方程

啮合域

在这一步我们离散化系统域成更小的领域,或元素。离散化让我们代表的几何域和近似的解决方案在每个元素更好地代表整个域的解决方案。元素的数量决定了问题的规模。

推导级方程式

在这一步中,我们假设的近似解微分方程在一个元素在指定点,或节点。在我们的示例中,解决方案确定的离散值的\φ(\ \)位移\ (x \)和\ (y \)方向(2 d分析)。未知的主字段的数量在一个节点被称为节点的自由度(自由度)。

现在控制微分方程应用于单个元素的域(图6)。在元素水平,控制方程的解由一个连续函数近似代替\φ(\ \)的分布在元素域\ (D e ^ \),表达的未知节点值\ \ phi_1 \, \ \ phi_2 \,和\ (\ phi_3 \) \φ(\ \)的解决方案。一个方程组的\ \ phi_1 \, \ \ phi_2 \,和\ (\ phi_3 \)可以制定的元素。

组装

我们获得解决方程结合解方程组的每个元素,以确保在每个节点连续性。在我们的示例中,元素刚度矩阵和力(\ (K_e \) \ (F_e \))创建一个全球刚度矩阵组装(\ (K \))和一个力矩阵(\ \)(F)在整个域

施加边界条件

我们在边界节点施加必要的边界条件和解决全球方程组。解决方案\φ(x, y)(\ \)问题成为一个分段近似,表示的节点的值\φ(\ \)。一个线性代数方程组装配过程的结果:\ (K \φ= F \)。

非常不寻常的实际工程问题有一个系统包含成千上万的方程。并行计算技术可以大大减少所需的时间组装矩阵和计算问题的解决方案的巨大规模。

解的方程

我们解决全球方程组\ (K \φ= F \)的位移\ (x \)和\ (y \)方向。

并行的问题

我们实现每一步在MATLAB有限元方法中定义^®函数。网格是通过使用mesh2d,基于MATLAB的网一代申请2 d几何图形。MATLAB中的分析工具显示,最耗时的操作属于刚度矩阵组装步骤。

这一步涉及到三个主要业务为每个元素:
1。计算单元刚度矩阵(\ (K_e \))级解方程组。\ \ (K_e \)的大小(N_e N_e \ \倍)
在那里,

2。地图的局部位置\ (K_e \)矩阵值位置的全局刚度矩阵。
3所示。全球刚度矩阵填充(\ (K \))使用的地图元素刚度矩阵的值。

表1。装配时间和总执行时间的比较不同的景深。

随着元素的数量增加,如此迭代的数量和大小的\ (K \)。图7和表1比较大会所花费的时间和系统的总执行时间在串行模式下不同的景深。显然,装配是最耗时的一部分,占据了99%以上的总执行时间当系统有很高的自由度。

刚度矩阵组装

幸运的是,最后组装\ (K \)矩阵是独立的元素的顺序选择内循环。我们可以评估几个元素的贡献全球刚度矩阵(\ (K \))同时通过MATLAB计算分布到多个工人。通常在一个串行执行装配操作为通过每个元素循环,哪些步骤,并确定其对全球刚度矩阵的贡献。我们只是把串行为循环一个平行为循环使用parfor在并行计算工具箱构建™(图8)。

在我们的示例中,全球刚度矩阵元素刚度矩阵的贡献的总和所有的迭代循环。的parfor减少结构使我们能够处理此类作业(通常的形式\ (r = f ({expr} \文本,r) \))自动。

演示性能通过使用并行的方法,问题是扩大从粗网格极其高雅。粗网格包含大约128元素共有150景深。我们细化网格,直到861122年它包含856800个元素与景深。我们细化网格,悬臂梁的自由端位移的聚合

并行的方法中,我们使用一个计算机集群中有一个头节点和5机器,每个使用以下配置:双英特尔^®至强^®1.6 GHz四核处理器(8个核每台机器总共40核心),13 GB RAM与Windows^®64位操作系统。每台机器跑8 MATLAB的工人总共40。测量串行执行时间,我们使用一个MATLAB工人头节点上运行。使用一个64位的操作系统使我们能够创建大型稀疏矩阵(861122 x 861122)没有运行到内存限制。在大多数有限元应用,合成K在本质上是稀疏的。

数字9和9 b比较总执行时间在秒增加自由度之间的串行(红色)和并行执行的(绿色)模式。表2总结了结果。

表2。代表数据用于数字9和9 b。

注意,与几个自由度,系统分配操作的成本远远高于这些操作的执行时间。因此,当系统只有250景深,串行执行实际上是速度比并行执行。

图9 b显示大约400景深,两条曲线相交的点,串行模式下执行速度比并行模式。我们看到后切换到并行模式的性能改进这一点。实际执行时间和交叉点取决于几个因素,包括MATLAB函数的执行速度,工人的机器的处理速度,网络速度,和工人数量,可用内存和系统负载。

总结

在本文中,我们演示了一个简单的并行有限元分析应用程序的方法。通过分析串行代码的性能,我们开始关注最计算密集型我们设置的一部分。使用简单的代码更改我们能够显著提高应用程序的性能,减少分析时间从60个小时800000自由度系统进行40-worker设置小于2小时。

¹工人从MATLAB分别MATLAB计算引擎运行会话。

²另外,天然橡胶Sumant Aluru和交流Cangellaris快速有限元建模的方法静电驱动MEMS。”国际期刊工程中的数值方法2009;77:1789 - 1808。