## 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 at 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 shows as:

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

Let say we have n dofs, and m of them are master dofs. If rigid elements do not connect two nodes with constrainted supports we can define a matrix Pf (m * n) which when is multiplied to the F will get a Fr vector which is force vector for __R__educed 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 whith 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 is three properties which RigidElement that should be setted:

- 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 s inside this.

### 3. AppliedLoadTypes

By defaults is empty, rigid element will be applied when structure is analysing with a LoadCase which have 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 is connected with a RigidElement member. Other structure also cosidered but instead of RigidElement, nodes within a floor diaphragm are connected whith elements with factor C of average stiffness of members. In below chart you can see in every C , difference between node displacements of two structures. As you can see second structure is converged to first structure in large C s.

## Example

for example want to calculate the Kr and Mr of below grid where every nodes in a roof are 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

v1, Dec 7 2014