# Finite Element Mass and Stiffness Matrices

I am attempting to use the finite element method to solve a vibration problem. I am following user21's answer from here which is very helpful. My problem is that I am confused by the use of the mass and damping matrix. Eventually I will need both. Since force is mass times acceleration I think the second derivative with respect to time should appear in the differential equation to represent the acceleration. The first derivative is damping. user21 uses the first derivative and treats it like mass. I attempt to rectify this below but I get a mass matrix that is empty and a damping matrix that is full. What is happening?

I start by defining the plane stress differential equations for stress

Needs["NDSolveFEM"];
(* plane stress equations *)
ClearAll[planeStress];
planeStress::usage =
"planeStress[u,v,t,x,y,Y,\[Nu]] Y is modulus of elasticity and \
\[Nu] is Poission ratio. ";
planeStress[u_, v_, t_, x_, y_,
Y_, \[Nu]_] := {Inactive[
Div][{{0, -((Y*\[Nu])/(1 - \[Nu]^2))}, {-(Y*(1 - \[Nu]))/(2*(1 - \
\[Nu]^2)), 0}}.Inactive[Grad][v[t, x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-(Y/(1 - \[Nu]^2)),
0}, {0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}}.Inactive[Grad][
u[t, x, y], {x, y}], {x, y}],
Inactive[Div][{{0, -(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2))}, \
u[t, x, y], {x, y}], {x, y}] +
Inactive[
Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)),
v[t, x, y], {x, y}], {x, y}]}


Next I make the mesh and set up the boundary conditions.

(* make mesh *)
\[CapitalOmega] = ImplicitRegion[True, {x, y}];
mesh = ToElementMesh[\[CapitalOmega], {{0, 5}, {0, 1}},
"MaxCellMeasure" -> 0.05];
(* Set up boundary conditions *)
bcs = DirichletCondition[{u[t, x, y] == 0, v[t, x, y] == 0}, x == 0];


Here I set up the PDE with the second derivative to give an acceleration term. In user21's formulation there is a first derivative.

(* Set up PDE *)
pde2D = {D[u[t, x, y], {t, 2}], D[v[t, x, y], {t, 2}]} +
planeStress[u, v, t, x, y, 10^3, 33/100] == {0, 0};


Now I start NDSolve and put in initial conditions for displacement and velocity

(* Start NDSolve *){state} =
NDSolveProcessEquations[{pde2D, bcs,
u[0, x, y] == 0, (D[u[t, x, y], t] /. t -> 0) == 0,
v[0, x, y] == 0, (D[v[t, x, y], t] /. t -> 0) == 0}, {u, v}, {t,
0, 1}, {x, y} \[Element] mesh,
Method -> {"PDEDiscretization" -> {"MethodOfLines",
"SpatialDiscretization" -> {"FiniteElement"}}}];


The next lines of code are standard to get a discretized model

(*Extract the finite element data:*)
femdata = state["FiniteElementData"];
initBCs = femdata["BoundaryConditionData"];
methodData = femdata["FEMMethodData"];
initCoeffs = femdata["PDECoefficientData"];
vd = methodData["VariableData"];
nr = ToNumericalRegion[mesh];
sd = NDSolveSolutionData[{"Space" -> nr}];
(*Discretize the PDE and the boundary conditions:*)
discretePDE = DiscretizePDE[initCoeffs, methodData, sd];
discreteBCs = DiscretizeBoundaryConditions[initBCs, methodData, sd];


Now I look at the various matrices

{load, stiffness, damping, mass} = discretePDE["SystemMatrices"]


This gives

Why is the stiffness and damping matrix filled? Should it not be the stiffness and mass?

The damping matrix looks like

MatrixPlot[damping]


What is happening?

Thanks

Top level functions like NDSolve, NDEigensystem will convert a second order (or higher) time derivative in a system of first order time derivative equations. The second order time derivative corresponds to the mass matrix and a first order time derivative corresponds to the damping matrix. This is done fully automatic and works for any time derivative.

Now, the reason that there is a mass matrix field in the data structure of the discretized PDE is to, for example, model Rayleigh damping. This needs to be done manually. There is more on that topic in the tutorial in the section A Swinging Beam—Transient Coupled PDEs.

Also, it is conceivable, to use a Verlet time integration scheme for PDEs of the form m.y''+d.y'+s.y==f.

Please have a look at the FEM Programming turorial which explains the structure in detail. Otherwise you may want to use NDEigensystem for you analysis which does this automatically.

• Thanks for your prompt comments. I am familar with converting a second order system in time to a first order system b
– Hugh
Nov 16, 2015 at 7:33
• (1) Thanks for your prompt comments.I am familiar with converting a second order system in time to a first order system by using intermediate variables.There are several approaches to this conversion particularly if you are including damping as well as mass in the starting differential equation.Does the damping matrix then always correspond to the coefficients of this converted first order system?What system of intermediate variables are used?
– Hugh
Nov 16, 2015 at 7:45
• (2) In your [example] (mathematica.stackexchange.com/a/57391/12558) the differential equation appears to be wrong (time should be second order) but is there a subtle trick I am missing? I am very familiar with the Rayleigh approximation but this is regarded as inadequate in careful work where damping must be represented more accurately.
– Hugh
Nov 16, 2015 at 7:47
• @Hugh, (1) yes, if you call ProcessEquations then the damping matrix will be the first order system. To get the intermediate variables look at the "VariableData` extracted from the "FEMMethodData" that should have the dependent variables. Concerning (2), note that the second order time derivative, the first order time derivative and the reaction term are all the same! The just end up in different matrices. Perhaps that's the confusion? Nov 17, 2015 at 7:28
• Thanks for the clarification. If I have a non-uniform distributed mass, a non-uniform distributed damping as well as a non-uniform set of material properties then how do I extract the matrix for each? Is this another question I should ask?
– Hugh
Nov 17, 2015 at 12:54