Understanding Finite-Element Modal Analysis

Question

I am teaching a basic course on computational physics and for the last part of the course I will introduce freshman physics undergraduates to finite-element modelling methods.

I am preparing a COMSOL model that solves for the eigenfrequencies and shape of the eigenmodes of the unperturbed 3D geometry of e.g. a tuning fork.

The problem is, I don't actually know what the FEM algorithm is doing: What matrices are set up, what goes into each element, what equations are solved for...?

The answer should give an account of what the FEM algorithm does to go from the meshed 3D geometry to the solved eigenfrequencies and the shape of the eigenmodes of the whole structure.

Answer 1

The biggest problem in answering to this question is that it's not as straightforward as you might think to see the physics in the equations until you see them fully developed, and even then you're doing your best to assign a meaning. As a graduate student in mechanical engineering, this was hard for me because I was used to using my physical intuition to guide me through the math forest.

For a FEM problem applied to modal analysis, you basically start with a strong form equation. For example, a free vibration of an elastic rod is given by: $$ \lambda \rho u + (Eu,_x),_x = 0,\ x \in (0, L) $$ Here, $\rho$ is the density, $E$ is the modulus of elasticity, and $u$ is the displacement. Let's apply a fixed boundary condition ($u(L) = 0$) and a free boundary condition (e.g. $-Eu,_x(0)=0$). You can transform it into the equivalent weak form by using a weighting function $w$: $$ \int_0^Lw,_xEu,_xdx\ -\lambda\int_0^Lw\rho u\, dx = 0 $$ This is where the math ends up losing it's physical interpretation. From here (and skipping a lot of rigor for brevity in this answer), we would basically restrict $w$ in a subset of the Sobolev space, apply the Galerkin approximation by discretizing the domain which creates the finite dimensional weighting function $w^h$ and solution $u^h$, make $u^h$ the sum of a Sobolev function $v^h$ and a given function $g^h$, and then define the functions $w^h$ and $v^h$ as a finite sum of linearly independent basis functions (i.e $v^h = \sum_{A=1}^N c_A N_A$ and $w^h = \sum_{B=1}^N d_A N_A$). By applying some linear algebra theorems and rearranging a lot, you eventually get the following equation: $$ (K-\lambda^h_kM)\psi_k = 0 $$ Here, $\psi_k$ is the kth eigenvector, which represents the displacement. It's hard to get an appreciation for the math that I've glossed over. The book The Finite Element Method has the full explanation for how this is developed (Chapter 1 goes over the basics of the method, and Section 7.3 develops the frequency analysis portion).

Just to clarify, we've simply called these matrices K and M, and the form of the equation bears a resemblance to a simple modal analysis for a beam, so we simply assign that as the physical interpretation. Really, the K and M matrices are highly dependent on the mesh and how the elements are connected. These large matrices are constructed from each element in the mesh.

The discretization is really important. Physically for this problem, the nodes represent points at which we're calculating the displacement, and the elements connected to that node restrict its movement under an applied force. You can also think of the mesh for this problem as having $n$ number of rods connected to each other. In that way, each rod piece only really affects those it's connected to. Each element is represented by a small matrix. The element matrices are calculated like this: $$ m_{ab}^e = \int_{\Omega^e} N_a \rho N_b\,d\Omega \\ k_{ab}^e = \int_{\Omega^e} N_{a,x} E N_{b,x}\,d\Omega $$ Here, $\Omega^e$ is the domain of the element itself. When the element matrices get assembled into the global matrix, some elements in the large matrix will have terms added from several smaller matrices. This is a representation of how the topology from the mesh affects the global matrix.

The $N_a$ and $N_b$ symbols are the basis functions that span the subset of the Sobolev space we chose. These are very important because different basis functions changes the solution. We can choose to use linear or parabolic functions, or even different types of functions, like Lagrange or even B-spline basis functions.

Again, I don't know how much FEM theory the OP has, and it's been really hard to try to sum up how all this works. This may not be sufficient, but hopefully it serves to highlight some of the theory behind the FE method.

Answer 2

Finite Element Analysis is a mathematical tool very extended among engineers. However, after more than a year researching on the topic of computer simulation, where FEA plays such an important role, I couldn't yet find a satisfactory explanation on how they really really work...

The main background of FEM is that of structural engineering in the 60s. Engineers are very practical people, so initially they devised a system which allowed them to set algebraic equations where the relations between different points or "nodes" of their structures could be set. These algebraic relations were further proven to be of many different types, not only structural, so also thermal relations could be formulated and many others: all one needed to analyse was a proper discretization of the space in the form of a mesh with nodes to relate to each other.

For our particular case of structural engineering the procedure can be summarized in its three main steps, so the main ideas come up in a graphic taste:

STEP 1: Discretization In this basic yet classic example, I have chosen to divide my sample structure in only 5 elements where nodes are clearly identifiable in the meetings between beams. There are four nodes:

(Degrees of Freedom)

Each node will have 6DOF (Degrees Of Freedom): three for linear displacement on each axis (X,Y,Z) and three for rotational around each axis (X,Y,Z), because we are working on 3D. Many examples available provide the more "simple" 2D situation, but in my opinion this only complicates things further.

Once the nodes are located and there is a network of how they relate to each other, we can consider we have a mesh. In our example, the correlation is depicted in the table: The elements serve to "link" nodes to each other. The table establishes a "topology" for the nodes.

STEP 2: Element characterization and shape functions In this step is where FEM formulation and literature get really really awkward and nasty. In fact, this is the core of everything and where FEM differs from other ways of solving PDEs.

In reality, and despite its mathematical complexity, what we are looking for is a way of characterizing the material properties of the element. For such purpose, the method requires that the behaviour of those links among nodes obeys some formula. This formula is the actual Shape Function. In fact, the shape function can be any mathematical formula that helps us to interpolate what happens wherever there are no points to define the mesh. This "ghost" entity that appears between nodes is in fact the Finite Element. In practical terms, as engineers we are more interested in the implementation of the FEM, not so much on its formulation, so what is important to understand is that for different shape functions we obtain different element matrices.

2DOF Beam element matrix 3DOF Beam Element matrix 3DOF Timoshenko Beam element matrix 6DOF Timoshenko Beam element matrix

Depending on the chosen formulation we have different degrees of interpolation and hence presumably higher or lower degrees of precision. Also depending on the chosen formulation we might have different ways of locating and relating our nodes to each other.

For our example we can choose any of the formulations provided in literature (above are the most common used in structural engineering). It is important to note the internal structure of these element matrices, which are symmetrical and clearly divided into parts each corresponding to the nodes that reside on the element's boundaries (2 nodes in the case of a beam - 4 quadrants in the matrix).

STEP 3: Matrix assembly and solution

Because the relations between nodes need to be accomplished all at the same time, we have to set all the equations in such a manner that they compose an algebraic system of equations. The matrix equation we want to solve (at least in statics) is as follows: \begin{equation} [F]= [Kg]·[u], \end{equation}

Where [F] is the vector of applied external forces, [Kg] is the system's global stiffness matrix and [u] is the vector of displacements resulting from the application of the forces. The size of [Kg], [F] and [u] is that of the number of DOFs times the number of nodes, being [Kg] squared and [F] and [u] unidimensional. In order to generate the Force vector, all we need to do is collect the applied forces (linear and moment forces) and sort them according to the index of the node they are applied to. For the global stiffness matrix, it is necessary a bit more laborious procedure by means of which we iterate throughout each element's particular stiffness matrix. Out of each one of those, we get only the part that corresponds to the position of the node we are storing in the matrix, and add it to the possible concurrent data that comes from other elements. Warning: before entering in the global stiffness matrix, we must convert local coordinates to global coordinates!

To get the displacement vector, it is needed to first enforce the constraints and then solve the resulting algebraic system of equations. To do this there are two classical approaches: Penalty Method and Lagrange Multipliers method, but we wont enter into this here...

Global Matrix Assembly

Afterwards, all we need to do is use any available algebraic equation solver (LU decomposition is one of the most extended), and obtain the solution of the system.

STEP 4: Computation of the eigenvalues of a matrix

Once you have a matrix assembled and representing your physical model, you can proceed to make other mathematical operations on it, like for example the computation of its eigenvalues. The eigenvalues have the "magical" property of, being the stiffness matrix a form of transformation of the displacement vector, represent that form in which the displacements are the largest. These displacements are directly linked with the physical properties of the structural system and its period of free vibration.

Answer 3

If you are familiar with the standard FEM analysis works, the idea of modal analysis is straightforward. In standard FEM analysis, you transfer the time-dependent elastic wave equation (ignoring damping for now) $$ \ddot{r}(t) + A r(t) = f(t) \tag{1}\label{1},$$ which is a mathematical model describing the behavior of the displacement $r(t)$ using the differential operator $A$ (which includes the material law; in the case of linear elasticity, $Ar = \mu \Delta r + (\mu+\lambda)\nabla(\nabla\cdot r)$, where $\lambda$ and $\mu$ are the Lamé parameters) into the system of linear equations $$ M \ddot{\vec{r}}(t) + K \vec r(t) = M\vec f(t),$$ where $K$ is the stiffness matrix and $M$ is the mass matrix, as you write in your question. The way this is done in finite element methods is by projecting the infinite-dimensional system \eqref{1} onto a finite-dimensional subspace (this is called Galerkin approach), which in finite element methods is chosen as the span of piecewise polynomials on each cell of the mesh. There are already a lot of answered questions about how this works on this site (e.g., this one), so I'm not going to go into more detail here; I'll just point out that in the specific case of piecewise linear functions, the finite-dimensional approximation to the solution of \eqref{1} is uniquely determined by its values at the nodes of the mesh, and this is exactly what $\vec r$ represents.

Similarly, the eigenmodes -- which are stationary states! -- satisfy the resonance equation $$ Ar = \lambda r,$$ where $\lambda$ is the eigenvalue. Proceeding as for the time-dependent case \eqref{1}, you arrive at the finite-dimensional (generalized) eigenvalue problem $$ K\vec r = \lambda M\vec r.$$ This is now a standard problem in numerical linear algebra, for which several methods exist (which, it should be pointed out, work by primarily computing eigenvectors, from which the eigenvalues can be obtained easily). In particular, since $K$ and $M$ are large and sparse, a Krylov method such as Arnoldi would be most efficient. This is exactly what COMSOL does (using the FORTRAN library ARPACK, which implements these methods).

(You are right that eigenvectors are not unique; usually -- and in particular in Krylov methods -- you want to have them normalized, to $r^TMr=1$ as Bill Greene notes, see page 55 of the COMSOL manual; this only leaves the choice of sign, which is arbitrary and often even random (due to the initialization of the method).)