# How to compute overall inertia properties from FE mass matrix?

How should I evaluate the overall mass and moments of inertia (polar and transverse) of a finite elements model, having its mass matrix?

Given that the global mass matrix is composed of symmetric mass submatrices (let's call these $$m_{i,j}$$, with $$i,j = 1,2,3,\ldots, N$$, $$N$$ = [number of nodes]) carrying information on inertia properties of each single node, like this:

$$\left[ \array{ M & & & & & \\ & M & & & \mathrm{symmetric} & \\ & & M & & &\\ & & & I_{11} & & \\ & & & -I_{21} & I_{22} & \\ & & & -I_{31} & -I_{32} & I_{33}\\ }\right]$$

I tried for example computing overall mass by summing mean($$m_{i,j_{1,1}}$$,$$m_{i,j_{2,2}}$$,$$m_{i,j_{3,3}}$$) for all $$i = j$$ (i.e. submatrices lying on the diagonal of the global matrix). This, however, yields wrong results.

• I assume that the submatrices in your global mass matrix above are not diagonal? Is there something special about the elements you want to consider? In a general 3D FE model, where the element mass matrices are calculated using a "consistent" formulation, the global mass matrix will have non-zero elements far off the diagonal; not at all like the version you show above. Feb 25 at 0:44
• Sorry if it took time for me to answer... Btw you're right, I am referring to a mass matrix of a (dynamically) reduced order model Mar 7 at 23:18
• OK, I supplied an answer that I think will do what you want. Mar 8 at 0:01

This note describes how to compute equivalent rigid body mass properties from a "mass matrix" produced by spatial discretization (finite element, finite difference, etc.) of a continuum. By equivalent, we mean that the kinetic energy computed for the equivalent rigid body will be equal to that for the mass matrix for any rigid body motion. We don't assume any particular structure for the mass matrix, e.g. it may be computed by a so-called consistent finite element formulation and therefore exhibit considerable coupling between degrees of freedom.

We assume a fixed Cartesian coordinate system whose origin is probably not at the center of mass. Often, the center of mass is something we want to compute.

The kinetic energy of the continuum is given by

$$E_c = \frac{1}{2} \int\limits_V \rho \dot u^T\cdot \dot u dv$$ where $$\rho$$ is the material density and $$u$$ is the vector of three displacements defined at every point in the domain of the continuum, $$V$$. The kinetic energy based on the discrete mass matrix is given by $$E_d = \frac{1}{2} \dot U^T M \dot U$$ where $$M$$ is the mass matrix and $$U$$ is the vector of displacements at each node in the model. This vector often consists of the three displacements at each node but may, in addition, include the three nodal rotations as degrees of freedom ($$\theta_x, \theta_y, \theta_z$$). For generality, we will consider this second case. Finally, the kinetic energy of the equivalent rigid body is given by $$E_r = \dot q^T \bar M \dot q$$ Where $$q$$ is the vector of six rigid body motions (three displacements and three rotations) relative to the reference coordinate system. $$\bar M$$ is the matrix we want to calculate. $$\bar M$$ is symmetric and for all practical cases the $$(1,1), (2,2), (3,3)$$ terms are all equal the total mass. The lower, right $$3\times 3$$ sub matrix is the inertia tensor

$$\left[ \begin{array}{ccc} I_{xx} I_{xy} I_{xz} \\ I_{xx} I_{xx} I_{xz} \\ I_{zx} I_{zy} I_{zz} \end{array} \right]$$

The following terms are of interest because they can be used to calculate the center of mass

\begin{align} \bar M(2,6) = \int\limits_V \rho x dv \\ \bar M(3,4) = \int\limits_V \rho y dv \\ \bar M(1,5) = \int\limits_V \rho z dv \end{align}

The displacement vector from the discrete model, $$U$$, approximates the displacements of the continuum, $$u$$. However, to calculate $$\bar M$$ we must restrict these displacements to only those representing rigid body motions. This can be seen as a change of basis, $$U=Qq$$ where $$Q$$ is a matrix of six basis vectors. The first three vectors are unit displacements in the $$x,y,z$$ directions, respectively and the 4th, 5th, and 6th vectors are the linearized rotations about the $$x,y,z$$ axes. (by linearized we mean that $$sin(\theta)\approx\theta$$ and $$cos(\theta)\approx1$$.) Constructing the first basis vector, for example, just requires setting the x-displacement at all node points to one and all other entries to zero. The second and third basis vectors are constructed analogously. Constructing the fourth basis vector requires setting all $$\theta_x$$ nodal degrees of freedom to one, $$v_i=-z_i$$ and $$w_i=y_i$$. Here $$v_i$$ and $$w_i$$ indicate the y-displacement and z-displacement at the ith node, respectively. The fifth and sixth basis vectors are constructed analogously.

With this basis change $$\bar M$$ is seen to be $$\bar M = Q^T M Q$$

From this $$\bar M$$ the mass, moment of inertia tensor, and center of mass can be directly identified as described above. And, of course various transformations can be applied to these terms as described in standard mechanics texts (e.g. calculate the inertia tensor about the center of mass).

Although this is not a direct answer to your question, what one should know here is how inertia moment calculations can be performed by FE in general.

For this purpose and in order to avoid trivial technicalities, I will assume that one has the mass matrix for a scalar FE model, i.e. a model with only 1 dof per node. This basic mass matrix is defined as

$$M_{ij}=\int_\Omega N_i N_j~d\Omega$$

where $$N_i$$, $$N_j$$ are shape functions of dofs $$i$$ and $$j$$.

All inertia related quantities are integrals which are rather trivial to discretize following standard FE theory.

For example, the total volume is defined as

$$V=\int_\Omega d\Omega$$

and using the partition of unity property, $$\sum_{i,j} N_i N_j = 1$$, the latter can be rewritten as $$V=\sum_{i,j} \int_{\Omega_e} N_i N_j~d\Omega =\begin{bmatrix}1&1&\dots&1\end{bmatrix}M \begin{bmatrix}1\\1\\\vdots\\1\end{bmatrix}$$

Knowing the $$x$$-coordinate $$x_i$$ of each dof $$i$$, the $$x$$-coordinate of the center of mass can be found as

$$x_c=\dfrac{1}{V}\int_\Omega x~d\Omega =\dfrac{1}{V}\sum_{i} \int_{\Omega} x_i N_i~d\Omega =\dfrac{1}{V}\sum_{i,j} \int_{\Omega} x_i N_i~N_j d\Omega =\dfrac{1}{V}\begin{bmatrix}x_1 & x_2 & \dots & x_N \end{bmatrix}M \begin{bmatrix}1\\1\\\vdots\\1\end{bmatrix}$$

which again is based on the partition of unity property $$\sum_{j} N_j = 1$$.

Assuming that the center of mass coordinates $$x_c$$,$$y_c$$,$$z_c$$ have been computed like this, they can be used for calculating moments of inertia. For example

$$I_{xz}=\int_\Omega (x-x_c)(z-z_c)~d\Omega =\dots =\begin{bmatrix}x_1-x_c & x_2-x_c & \dots & x_N-x_c \end{bmatrix}M \begin{bmatrix}z_1-z_c\\z_2-z_c\\\vdots\\z_N-z_c\end{bmatrix}$$

In general, you just need to rewrite the definition integral of the requested quantity in terms of the known partial integrals available in your mass matrix. Here I have shown the simplest case, in order to make clear the concept, but once you understand this methodology you can adapt it to any system.