VB.NET 有限元分析
4.93/5 (23投票s)
用于在 VB.NET 中完整实现有限元分析的强大框架
引言
关于有限元方法的说明在许多书籍中都有介绍。免费的初级代码在书籍和互联网上都很容易获得。然而，提供对有限元过程实现进行全面阐述的程序却很少见。此外，完全用 VB.NET 编写的有限元程序也不容易找到。
本文档试图为执行二维平面应力/平面应变问题的有限元分析提供一个全面的框架。所选用的单元虽然简单，但足以表达构建有限元程序的概念。
此外，通过利用方程的带状特性，可以高效地求解有限元方程组。虽然一些书籍确实提供了带状求解器的代码，但它们可能效率不高，或者编写的风格不甚清晰。本文还将介绍一个强大的并行带状求解器来填补这一空白。
预计该框架将为程序员编写自己的有限元分析程序提供一个良好的起点。
背景
有限元方法是一门非常成熟的科学。本文无意解释该过程背后的数学原理。建议读者在尝试编程之前，充分理解该方法。虽然会解释一些基本概念，
有限元方法
让我们快速回顾一下有限元方法的基本原理。请记住，除非您彻底理解并欣赏其基本思想，否则几乎不可能对有限元方法进行编程。
大多数物理现象可以用偏微分方程来表示，通常是高阶的。但求解这些微分方程可能无法实现或非常繁琐。在这种情况下（通常是这样），我们诉诸数值方法来求解。有限元方法就是这样一种数值方法。有限元方法 (FEM) 将偏微分方程表示为一组线性代数方程。应该注意的是，在将偏微分方程转换为一组线性代数方程时，可能会丢失显著的精度。因此，有限元分析的结果存在固有的误差。此外，计算过程也会给解决方案带来误差。因此，在接受任何有限元分析的结果作为问题正确解决方案之前，必须仔细研究。
有限元解决方案的一个重要部分涉及单元的开发。核心来说，单元表征了溶液在整个域上的变化。溶液的质量在很大程度上取决于所使用的单元的质量。
有限元分析步骤
一般来说，任何有限元分析都需要执行以下步骤。
- 离散化连续介质
- 生成单元矩阵：刚度矩阵；载荷向量
- 组装单元矩阵：形成全局刚度矩阵
- 施加边界条件
- 求解方程并计算变形
- 计算应变
- 计算应力
连续介质的离散化涉及将给定的问题域划分为若干子域，称为“单元”。每个单元都由称为“节点”的虚点界定和定义。例如，考虑下图所示的梁。
为了用有限元方法分析梁，我们必须细分该域（离散化）为若干有限单元。单元的类型是分析者必须根据其工程经验和判断选择的。如果我们用四边形单元来离散化域，可能会是这样。
单元由节点定义。在域中，每个单元可以单独视为
离散化或“网格划分”，这是更常见的说法，它本身是一门庞大的科学，在本篇文章的讨论范围之外。
每个单元假设单元域内的位移场。基于此假设，可以生成单元的刚度矩阵。这里假设使用一个具有单元域内线性变形的 3 节点三角形单元。该单元具有线性的形状函数，在每个节点处具有水平和垂直平动。节点 1 处的形状函数和每个节点的自由度可以表示为
该单元非常基础和简单，但足以提供基本理解。该单元每个节点有 2 个自由度，因此产生一个 6x6 的刚度矩阵。由于变形是线性的，因此应变和应力都是恒定的。因此，该单元通常被称为恒应变三角形 (CST)。刚度矩阵可以通过最小势能法、Rayleigh-Ritz 法或其他变分法推导出来。请读者参考更多关于该主题的基础文献。为了完整起见，此处列出了推导单元刚度矩阵的关键方程。
每个节点有两个自由度——u 和 w，或者通常称为 u 和 v。形状函数近似了场在单元域上的变化。因此，在单元的任何点处，位移可以表示为
其中，[N] 是形状函数矩阵，形状函数由以下公式给出：
其中，A 是三角形单元的面积，可以通过雅可比矩阵行列式的一半来确定，或者表示为
生成刚度矩阵的常用过程涉及计算雅可比矩阵 **[J]**、应变-位移连接矩阵“**[B]**”以及最终计算刚度矩阵 **[Ke]**，如下所示：
其中，根据常用表示法，**[B]** 是应变-位移连接矩阵，**[D]** 是本构矩阵。
如果 \(x_i\) 和 \(y_i\) 是给定单元节点 \(i\) 的坐标，则 \(x_{ij}\) 定义为 \(x_{ij} = x_i - x_j\)。应变的近似和随后的应变-位移矩阵 **[B]** 的计算可以如下进行：
二维膜问题中的本构矩阵取决于平面应力或平面应变条件的假设。对于平面应力条件，[D] 矩阵如下所示：
对于平面应变条件：
其中，\(E\) 是材料的弹性模量（各向同性），\(\nu\) 是泊松比。雅可比矩阵 \(\left[J\right]\) 如下：
三角形的面积由 \(A = \frac{1}{2} |J|\) 给出。
[B] 矩阵是基于二维坐标（x 和 y）的。因此，乘积 \(B^T D B \) 给出了体积积分的面积部分。刚度矩阵的最终方程为：
其中，t 是单元的厚度。
上述生成单元刚度矩阵的程序可以在代码中表示为：
''' <summary>
''' Returns the determinant of the Jacobian
''' </summary>
''' <returns></returns>
Private Function getDetJ() As Double
Return (x13 * y23 - x23 * y13)
End Function
''' <summary>
''' Returns the [B] matrix
''' </summary>
''' <returns></returns>
Public Function getB() As Double(,)
Dim B(2, 5) As Double 'initialize a 3x6 array
Dim detJ As Double = getDetJ()
B(0, 0) = y23 / detJ
B(1, 0) = 0.0
B(2, 0) = x32 / detJ
B(0, 1) = 0.0
B(1, 1) = x32 / detJ
B(2, 1) = y23 / detJ
B(0, 2) = y31 / detJ
B(1, 2) = 0.0
B(2, 2) = x13 / detJ
B(0, 3) = 0.0
B(1, 3) = x13 / detJ
B(2, 3) = y31 / detJ
B(0, 4) = y12 / detJ
B(1, 4) = 0.0
B(2, 4) = x21 / detJ
B(0, 5) = 0.0
B(1, 5) = x21 / detJ
B(2, 5) = y12 / detJ
Return B
End Function
Private Function getAe() As Double
Return getDetJ() * 0.5
End Function
Public Function getD() As Double(,)
Dim D(2, 2) As Double
Dim c As Double
If ProblemType = ProblemTypes.PlaneStress Then
c = ElasticityModulus / (1.0 - (PoissonRatio * PoissonRatio))
D(0, 0) = 1.0 * c
D(1, 0) = PoissonRatio * c
D(2, 0) = 0.0
D(0, 1) = PoissonRatio * c
D(1, 1) = 1.0 * c
D(2, 1) = 0.0
D(0, 2) = 0.0
D(1, 2) = 0.0
D(2, 2) = (1 - PoissonRatio) * 0.5 * c
Else
'plane strain
c = ElasticityModulus / ((1 + PoissonRatio) * (1.0 - 2 * PoissonRatio))
D(0, 0) = (1.0 - PoissonRatio) * c
D(1, 0) = PoissonRatio * c
D(2, 0) = 0.0
D(0, 1) = PoissonRatio * c
D(1, 1) = (1 - PoissonRatio) * c
D(2, 1) = 0.0
D(0, 2) = 0.0
D(1, 2) = 0.0
D(2, 2) = (0.5 - PoissonRatio) * c
End If
Return D
End Function
Public Function getKe() As Double(,)
Dim B(,) As Double = getB()
Dim BT(,) As Double = getBT(B)
Dim D(,) As Double = getD()
Dim teAe As Double = Thickness * getAe()
'Compute Ke now
Dim BTD(,) As Double = MultiplyMatrices(BT, D)
Dim Ke(,) As Double = MultiplyMatrices(BTD, B)
'multiply ke by teAe
Dim i, j As Integer
For i = 0 To 5
For j = 0 To 5
ke(i, j) = ke(i, j) * teAe
Next
Next
Return ke
End Function
Private Function getBT(ByRef B(,) As Double) As Double(,)
'returns the transpose of [B]
Dim BT(5, 2) As Double
Dim i, j As Integer
For i = 0 To 2
For j = 0 To 5
BT(j, i) = B(i, j)
Next
Next
Return BT
End Function
''' <summary>
''' Multiplies matrix a by b and returns the product
''' </summary>
''' <param name="a"></param>
''' <param name="b"></param>
''' <returns></returns>
Private Function MultiplyMatrices(ByRef a(,) As Double, ByRef b(,) As Double) As Double(,)
Dim aCols As Integer = a.GetLength(1)
Dim bCols As Integer = b.GetLength(1)
Dim aRows As Integer = a.GetLength(0)
Dim ab(aRows - 1, bCols - 1) As Double
Dim sum As Double
For i As Integer = 0 To aRows - 1
For j As Integer = 0 To bCols - 1
sum = 0.0
For k As Integer = 0 To aCols - 1
sum += a(i, k) * b(k, j)
Next
ab(i, j) += sum
Next
Next
Return ab
End Function
组装方程
单元方程给出了单元的属性，而不是被分析的体的属性。要分析该体，必须将所有有限元连接在一起。这个过程在数学上称为组装方程。每个单元刚度矩阵都被组装成一个全局矩阵。单元矩阵是一个大小为（单元 DOFs x 单元 DOFs）的方阵。整个体的刚度矩阵大小为（总 DOFs x 总 DOFs）。通常，组装体的总自由度可能在几十万的量级。我们通常会尝试以一种方式进行网格划分，使总自由度尽可能少（几千），但对于复杂的结构，出现几百万的自由度是很常见的。那么，一个中等规模的 10,000 DOF 结构的内存需求将是：
矩阵中的条目数 = \(10,000 \times 10,000 = 10^8\)。
每个条目如果保存为双精度数，则需要 8 字节。因此，整个矩阵将需要 **0.8 GB**。操作如此大的矩阵将非常耗时且难以进行。对于中等规模的结构，我们可能会得到一个无法容纳在系统内存中的全局刚度矩阵。然而，刚度矩阵是稀疏的，包含 **大量** 的零。这些零可以避免进行计算。这允许我们使用智能方案来存储刚度矩阵。刚度矩阵通常以以下格式存储：
- 带状矩阵
- 轮廓线存储
该框架将全局刚度矩阵存储为带状矩阵形式。方程求解器也必须相应修改，以处理所采用的存储方案。本文实现并提供了基于带状矩阵的 Gauss 消元求解器。该求解器在因子分解块上实现并行处理，使其快速而健壮。
单元的组装涉及根据自由度将每个单元矩阵放置在全局矩阵的适当位置。通过一个例子可以最好地理解这个过程。
考虑下图所示的悬臂梁示意图。
一个好的网格划分是将该梁细分为大量单元。然而，为了理解这个过程，我们将该梁细分为一个粗略的网格，如下所示。
每个节点处的自由度由节点编号绑定。因此，每个节点处的自由度可以表示为：
为每个单元形成的单元刚度矩阵将因此被放置在全局刚度矩阵中。下面显示了一个例子。
这就产生了一个组装好的全局刚度矩阵，如下所示：
这里描绘的算法实现如下：
Private Sub AssembleKg(ByRef Ke(,) As Double, ByRef Kg(,) As Double, ElementNo As Integer)
Dim i, j As Integer
Dim dofs() As Integer = {getDOFx(Elements(ElementNo).Node1 - 1),
getDOFy(Elements(ElementNo).Node1 - 1),
getDOFx(Elements(ElementNo).Node2 - 1),
getDOFy(Elements(ElementNo).Node2 - 1),
getDOFx(Elements(ElementNo).Node3 - 1),
getDOFy(Elements(ElementNo).Node3 - 1)}
'Place the upper triangle of the elemental stiffness matrix in the global
'matrix in proper position
Dim dofi, dofj As Integer
For i = 0 To 5 'each dof of the ke
dofi = dofs(i)
For j = 0 To 5
dofj = dofs(j) - dofi
If dofj >= 0 Then
Kg(dofi, dofj) = Kg(dofi, dofj) + Ke(i, j)
End If
Next
Next
End Sub
Private Function getDOFx(NodeNo As Integer) As Integer
Dim nDofsPerNode As Integer = 2
Return (NodeNo) * nDofsPerNode
End Function
Private Function getDOFy(NodeNo As Integer) As Integer
Dim nDofsPerNode As Integer = 2
Return NodeNo * nDofsPerNode + 1
End Function
Private Function getHalfBandWidth() As Integer
Dim i As Integer
Dim n(2) As Integer
Dim MaxDiff As Integer = 0
Dim diff As Integer
For i = 0 To NElements - 1
n(0) = Elements(i).Node1
n(1) = Elements(i).Node2
n(2) = Elements(i).Node3
diff = n.Max - n.Min
If MaxDiff < diff Then MaxDiff = diff
Next
'now we have maxdiff
'half band width is maxdiff * 2. 2 because there are 2 dofs per node
Dim hbw As Integer = (MaxDiff + 1) * 2
If hbw > NNodes * 2 Then hbw = NNodes * 2
Return hbw
End Function
请注意，如果我们交换节点 1 和 2，编号将更有效率，结果的带宽将更紧凑。上面显示的矩阵的轮廓线存储如下：
轮廓线存储方案显然更健壮，尽管如此，本文采用了半带宽存储方案。深紫色方块显示非零项。因此，避免零项可以让我们只存储如图所示的带状区域。
此外，线性问题的刚度矩阵本质上是对称的。这允许我们只存储带宽的一半，但仍然得到正确的解。这就需要采用本文框架中采用的半带宽存储方案。
刚度矩阵描述了结构/体在空间中的浮动状态。此时的刚度矩阵是奇异的，意味着如果结构受到轻微扰动，它将表现出无限的刚体变形。为了进行有意义的分析，我们必须通过施加边界条件来约束该体。边界条件可以是载荷和支撑的形式。虽然载荷可以有很多类型，但本文仅以节点力为例。此外，支撑也可以有很多类型。只考虑 x 和 y 方向的约束。
通过惩罚方法施加支撑。这包括将与受约束自由度对应的对角线项乘以一个大数。关于惩罚方法背后深入数学的讨论被认为超出了本文的范围。
一旦施加了边界条件，下一步就是使用 Gauss 消元法求解方程并计算每个自由度处的变形。该框架在带状矩阵上执行 Gauss 消元。在计算变形后，可以提取单元变形，并通过上述公式 \( \{ \epsilon \} = [ B ] \{ d \}\) 计算应变。然后根据胡克定律计算单元应力。
绘制单元应变或应力需要定义一个颜色映射。颜色映射定义了从最小值到最大值的颜色变化。例如，如果应力范围是 -20 到 350，颜色映射将在 -20 和 350 之间分配指定的颜色。完成颜色分配后，通过线性插值可以轻松确定所需值对应的颜色的索引。这里，颜色映射被定义为在红-黄-绿-青-蓝之间变化。创建颜色映射的代码如下所示：
Private Sub setColors()
'Dim mainColors() As Color = {Color.Red, Color.Yellow, Color.LightGreen,
' Color.Cyan, Color.Blue}
Dim mainColors() As Color = {Color.Blue, Color.Cyan, Color.LightGreen,
Color.Yellow, Color.Red}
Dim nColors As Integer = mainColors.Length
'this is from top to bottom.
Dim percentStep As Double = 100 / (nColors - 1)
Dim Indices(nColors - 1) As Integer 'we will have these many main indices
Dim i, j As Integer
For i = 0 To nColors - 1
Indices(i) = CInt(i * percentStep / 100 * nSegs)
Next
Indices(0) = 0
Indices(Indices.Length - 1) = nSegs
ReDim Colors(nSegs)
Dim c() As Color
For i = 0 To Indices.Length - 2
c = getSteppedColors(mainColors(i + 1), mainColors(i), Indices(i + 1) - Indices(i))
For j = 0 To c.Length - 1
Colors(Indices(i) + j) = c(j)
Next
Next
End Sub
一旦颜色在数组中设置好，就可以简单地通过以下方式检索所需的颜色：
Public Function getColor(Value As Double) As Color
'value exists somewhere between min and max..
'we find color between min color and max color corresponding to value and return
Dim rangeVal As Double = Max - Min
Dim dv As Double = Value - Min
Dim ratio = dv / rangeVal
Dim colIndex As Integer = CInt(ratio * nSegs)
Return Colors(colIndex)
End Function
使用代码
有限元框架是用 VB.Net 实现的，名为 btFEM。实现的主要类包括：
- CST - 恒应变三角形的实现。该类返回 **[B]**、**[D]** 和 **[Ke]**。
- btGauss - 用于求解方程 \( [K]{d} = {r}\) 的带状 Gauss 求解器的并行和串行实现。
- Element - 单元的实现。这与 CST 不同。Element 类包含上层信息，例如每个单元的边界节点编号等。
- Node - 包含每个节点坐标的定义。
- PointLoad - 包含指定节点在 x 和 y 方向上的载荷定义。
- ColorMap - 实现颜色条以显示体上的应力和应变变化。
启动文件：frmbtFem 实现打开模型的图形显示。请注意，必须将模型的网格提供给该程序。左键单击放大，右键单击缩小模型视图。滚动条可以平移模型，屏幕右下角的按钮可以将缩放重置为 1，平移值重置为 0。打开的有限元模型如下所示：
要分析结构，请点击菜单“模型”->“分析”。分析后，程序将显示变形后的形状，如下所示：
可以查看各种结果，例如 x 方向的应力、y 方向的应力和剪应力，或者 x 方向的应变、y 方向的应变和剪应变。例如，x 方向应力的图如下所示：
通过点击菜单“模型”->“结果”可以在文本框中查看变形、应力和应变。通过将变形和应力与标准问题（如自由端带有集中载荷的悬臂梁、简支梁）进行比较，该框架已经得到了验证。相同的输入文件已包含在项目文件中。下图显示了一个带孔板的典型验证运行，显示了应力 \(\sigma_x\) 的变化。
关注点
有限元分析是一项具有挑战性且计算量大的仿真。核心挑战在于高效计算单元矩阵和求解方程。Gauss 消元法是一种通用的方程求解技术，可以轻松高效地修改以处理带状矩阵。
本文描述了有限元方法的编程。提供完整的可运行程序有望增强读者的理解。健壮的带状组装器和求解器已作为独立类包含在内，可以轻松集成到任何代码中。
历史
文章版本 2。