Finite Element Mass and Stiffness Matrices

Question

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["NDSolve`FEM`"];
(* 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))}, \
{-((Y*\[Nu])/(1 - \[Nu]^2)), 0}}.Inactive[Grad][
      u[t, x, y], {x, y}], {x, y}] + 
   Inactive[
     Div][{{-(Y*(1 - \[Nu]))/(2*(1 - \[Nu]^2)), 
       0}, {0, -(Y/(1 - \[Nu]^2))}}.Inactive[Grad][
      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} = 
  NDSolve`ProcessEquations[{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 = NDSolve`SolutionData[{"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

Answer 1

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.