RigidElements in BriefFiniteElement.NET

Introduction

This article is about how to consider rigid bodies (non deformable elements, like rigid diaphragm) in the BriefFiniteElement.NET library which is an opensource library for linear finite element analysis in .NET platform.

Download BriefFiniteElement.NET Code

Please look here. After loading the "SOURCE CODE" tab, you have to select the 'RigidElement' as branch (by default 'master' branch is selected).

Background

In structural analysis, sometimes there is need to consider undeformable objects (element with infinite stiffness) in the FEM model. Examples are rigid diaphragms in earthquakes or rigid links for example for connecting a 1d element to a 2d element. Rigid elements have a very useful feature which is reducing the stiffness, mass and damp matrix dimensions and can significantly reduce the amount of calculations for analyzing the structure.

Formulation

*Note: this section is about formulation of rigid body elements. For using rigid element, you do not need to know anything about the formulation and this is just showing the formulation approach which library is using for handing rigid elements.

Since the models are linear, the master slave method is used for considering the rigid elements. You can find a detail about how it is with this paper. In linear analysis displacement vector, force vector and stiffness matrix should be determined and after some calculation, we should solve a linear equation system and displacements will be found.

Because of reducing the count of degree of freedoms (DoF) of structure, rigid element will reduce the stiffness, mass and damp matrix dimensions. For example, consider a problem with no settlement which can be shown as:

F = K . D   => D = (K^-1) * F

Let's say we have n dofs, and m of them are master dofs. If rigid elements do not connect two nodes with constrained supports, we can define a matrix  Pf (m * n) which when multiplied to the F will get a Fr vector which is force vector for Reduced structure (after applying the rigid elements to reduce DoFs). Also can define a Pd (n * m) matrix, which when is multiplied to Dr (Dr is displacement vector for reduced structure) will give the D which is displacement vector for original structure as below:

Fr = Pf * F

D = Pd * Dr

Can combine these two equations with first one as:

F = K . Pd . Dr  (pre multiply both sides with Pf)=> Pf . F = Fr = Pf . K . Pd . Dr

Taking Kr = Pf . K . Pd

Fr = Kr * Dr

This way, we can convert the problem into a reduced problem. Same can be applied for dynamic analysis:

M . Ẍ + C . Ẋ + K . X = F(t)

Taking:

X = Pd * Xr => Ẋ = Pd * Ẋr => Ẍ = Pd * Ẍr

Then:

M . Pd * Ẍr + C . Pd * Ẋr + K . Pd * Xr = F(t)

Premultiply which Pf:

Pf . M . Pd * Ẍr + Pf . C . Pd * Ẋr + Pf . K . Pd * Xr = Pf . F = Fr

Fr = Mr . Ẍr + Cr . Ẋr + Kr . Xr

where:

  Mr = Pf . M . Pd

  Cr = Pf . C . Pd

  Kr = Pf . K . Pd

This is the way that library is doing the calculation.

RigidElement Class

There is a RigidElement class, which is a rigid element! Using it is very simple, after creating a new RigidElement, should fill its Nodes property with nodes that rigid element connects together. Then define situation that rigid element should be applied. For example, in earthquake force that a single force is applying to the mass center of roof, the rigid diaphragm should be considered (if not, the elements which are connected to the mass center node will carry all lateral force of roof) and in for example live or dead loads rigid diaphragm should not be applied (because if apply, then all elements within rigid element will not have deformation and their internal force will be zero). For defining the situation, there are three properties which RigidElement should set:

• bool RigidElement.UseForAllLoads
• LoadCaseCollection AppliedLoadCases
• LoadTypeCollection AppliedLoadTypes

1. RigidElement.UseForAllLoads

It is false by default, if set to true then rigid element will be applied in every situation and all loads.

2. AppliedLoadCases

By defaults is empty, rigid element will be applied when structure is analysing with LoadCase inside this.

3. AppliedLoadTypes

By defaults is empty, rigid element will be applied when structure is analysing with a LoadCase which has a LoadCase, LoadType which is present inside this:

Validation

For validating the RigidElement class, a structure with randomized loading is used. Nodes in every floor are connected with a RigidElement member. Other structure also considered but instead of RigidElement, nodes within a floor diaphragm are connected with elements with factor C of average stiffness of members. In the below chart, you can see in every C, difference between node displacements of two structures. As you can see, the second structure is converged to the first structure in large C s.

Example

For example, I want to calculate the Kr and Mr of the below grid where every node in a roof is in a single rigid element, also all nodes of first floor are fixed to the ground. so will have 3 roofs = 3 * 6 = 18 DoF.

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)

Points of Interest

Anyone is welcome to fork the library and develop it...

History

• 7^th December, 2014: v1