RigidElements 简介 FiniteElement.NET
考虑刚体
引言
本文介绍了如何在 .NET 平台上用于线性有限元分析的开源库 BriefFiniteElement.NET 中考虑刚体（不可变形单元，如刚性楼板）的方法。
下载 BriefFiniteElement.NET 代码
请点击此处。加载“源代码”标签页后，您需要选择“RigidElement”作为分支（默认选择“master”分支）。
背景
在结构分析中，有时需要在有限元模型中考虑不可变形的物体（刚度无限的单元）。例如，地震中的刚性楼板或用于连接一维单元到二维单元的刚性连接件。刚性单元有一个非常有用的特性，那就是减小刚度、质量和阻尼矩阵的尺寸，从而显著减少结构分析的计算量。
公式
*注意：本节内容是关于刚体单元的公式推导。使用刚性单元时，您无需了解任何关于公式推导的知识，这里仅展示了库在处理刚性单元时所采用的方法。
由于模型是线性的，因此使用主从方法来处理刚性单元。您可以参考这篇论文来了解其详细信息。在线性分析中，需要确定位移向量、力向量和刚度矩阵，经过一些计算后，需要求解一个线性方程组来得到位移。
由于减少了结构的自由度（DoF）数量，刚性单元会减小刚度、质量和阻尼矩阵的尺寸。例如，考虑一个无沉降的问题，可以表示为：
F = K . D => D = (K^-1) * F
假设我们有 n 个自由度，其中 m 个是主自由度。如果刚性单元不连接两个具有约束支撑的节点，我们可以定义一个 Pf (m * n) 矩阵，该矩阵乘以 F 得到 Fr 向量，Fr 向量是减缩结构的（在应用刚性单元以减少自由度后）力向量。我们还可以定义一个 Pd (n * m) 矩阵，该矩阵乘以 Dr（Dr 是减缩结构的位移向量）得到 D（原始结构的位移向量），如下所示：
Fr = Pf * F
D = Pd * Dr
我们可以将这两个方程与第一个方程结合起来，得到：
F = K . Pd . Dr (两边同时乘以 Pf)=> Pf . F = Fr = Pf . K . Pd . Dr
令 Kr = Pf . K . Pd
Fr = Kr * Dr
这样，我们就可以将问题转化为一个减缩后的问题。动态分析也可以采用相同的方法：
M . Ẍ + C . Ẋ + K . X = F(t)
令
X = Pd * Xr => Ẋ = Pd * Ẋr => Ẍ = Pd * Ẍr
那么
M . Pd * Ẍr + C . Pd * Ẋr + K . Pd * Xr = F(t)
两边同时乘以 Pf
Pf . M . Pd * Ẍr + Pf . C . Pd * Ẋr + Pf . K . Pd * Xr = Pf . F = Fr
Fr = Mr . Ẍr + Cr . Ẋr + Kr . Xr
其中
Mr = Pf . M . Pd
Cr = Pf . C . Pd
Kr = Pf . K . Pd
这就是库进行计算的方式。
RigidElement 类
有一个 RigidElement 类，它就是一个刚性单元！使用起来非常简单，创建新的 RigidElement 后，需要用刚性单元连接的节点填充其 Nodes 属性。然后定义刚性单元应该应用的情况。例如，在地震力作用下，当单个力施加到屋顶的质心时，应考虑刚性楼板（否则，连接到质心节点的单元将承担屋顶的所有侧向力）；而在活荷载或恒荷载作用下，则不应应用刚性楼板（因为如果应用，则刚性单元内的所有单元将不会发生变形，它们的内力将为零）。为了定义情况，RigidElement 有三个属性需要设置：
bool RigidElement.UseForAllLoadsLoadCaseCollection AppliedLoadCasesLoadTypeCollection AppliedLoadTypes
1. RigidElement.UseForAllLoads
默认值为 false。如果设置为 true，则刚性单元将在所有情况和所有荷载下应用。
2. AppliedLoadCases
默认为空。当结构使用此集合中的 LoadCase 进行分析时，将应用刚性单元。
3. AppliedLoadTypes
默认为空。当结构使用包含此集合中的 LoadType 的 LoadCase 进行分析时，将应用刚性单元。
验证
为了验证 RigidElement 类，使用了具有随机加载的结构。每层楼的节点都用一个 RigidElement 成员连接。也考虑了其他结构，但不是用 RigidElement，而是用平均构件刚度因子 C 的单元连接楼层内的节点。在下面的图表中，您可以看到 C 值不同时两个结构节点位移的差异。正如您所看到的，在 C 值较大时，第二个结构收敛到第一个结构。
示例
例如，我想计算下图中网格的 Kr 和 Mr，其中屋顶上的每个节点都在一个刚性单元内，并且所有一楼节点都固定到地面。因此，将有 3 个屋顶 = 3 * 6 = 18 个自由度。
var nm = 4;
var str = StructureGenerator.Generate3DGrid(nm, nm, nm); //creating a structure
#region Create rigid elements
var roofs = str.Nodes.GroupBy(i => i.Location.Z).Skip(1).ToList();
foreach (var roof in roofs)
{
var nodes = roof.ToList();
var elm = new RigidElement(nodes.ToArray());
elm.UseForAllLoads = true;
str.RigidElements.Add(elm);
}
#endregion
var map = DofMappingManager.Create(str, LoadCase.DefaultLoadCase);//creating a map
var kt = MatrixAssemblerUtil.AssembleFullStiffnessMatrix(str);//assembling Kt
var mt = MatrixAssemblerUtil.AssembleFullMassMatrix(str); //assembling Mt
var Pf = PermutationGenerator.GetForcePermute(str, map); //calculating Pf matrix
var Pd = PermutationGenerator.GetDisplacementPermute(str, map);//Calculating Pd matrix
var kr = Pf * kt * Pd;//calculating Kr
var mr = Pf * mt * Pd;//calculating Mr
//Note: Kr contains both fixed and released parts (4 zones of Kff,Kfs,Ksf,Kss).
//Using CalcUtil.GetReucedZoneDevidedMatrix zones are separated automatically
var krDevided = CalcUtil.GetReucedZoneDevidedMatrix(kr, map);//dividing zones of Kr
var mrDevided = CalcUtil.GetReucedZoneDevidedMatrix(mr, map);//dividing zones of Mr
var krff = krDevided.ReleasedReleasedPart.ToDenseMatrix(); //free-free part of Kr,
//size: 18x18 (since there are 3 master nodes with each node 6 free DoF)
var mrff = mrDevided.ReleasedReleasedPart.ToDenseMatrix(); //free-free part of Mr,
// size: 18x18 (same as above)
关注点
欢迎任何人 fork 该库并进行开发...
历史
- 2014 年 12 月 7 日：v1