Generating 4-degree of freedom beam stiffness matrix
by Nasser M. Abbasi  oct 26,2009

This diagram below is from the book "Concepts and applications of finite element methods" 4th edition. The goal is to find the stiffness matrix for this beam. It has 2 degrees of freedom per node. I will start by assuming the deflection w(x) to be a approximated by a cubic polynomial as a function of x, and then use minimum potential energy approach to find the stiffness matrix.

Assume that the deflection is approximated by w=(Subscript[a, 1])+(Subscript[a, 2]) x+(Subscript[a, 3]) (x^2)+(Subscript[a, 4]) (x^3)and from this we need to obtain the shape function N so that we can write w=N (({{(Subscript[v, 1])}, {(Subscript[θ, 1])
(Subscript[v, 2])
(Subscript[θ, 2])}})) hence we see that N must be a 1 by 4 row vector. We call (Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[a, 4]) the generalized coordinates, and (Subscript[v, 1]),(Subscript[v, 2]),(Subscript[θ, 1,])(Subscript[θ, 2]) as the nodal degrees of freedom.

w=(Subscript[a, 1])+(Subscript[a, 2]) x+(Subscript[a, 3]) (x^2)+(Subscript[a, 4]) (x^3);

Now find an expression for each of the above 4 nodal unknowns in terms of the deflection. Notice that the angle θ=(dw/dx) as it is being approximated for small deflection (as in the angle of an arc (in radian) is the arc length over radius)

eq1=(Subscript[v, 1])==w/.x->0;
eq2=(Subscript[v, 2])==w/.x->L;
eq3=(Subscript[θ, 1])==D[w,x]/.x->0;
eq4=(Subscript[θ, 2])==D[w,x]/.x->L;

Now to solve for the generalized d.o.f (Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[a, 4])

generalDegreOfFreedomVars={(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[a, 4])};
First@Solve[{eq1,eq2,eq3,eq4},generalDegreOfFreedomVars];
wNew=w/.%;
Print["w(x)="<>ToString[wNew,FormatType->TraditionalForm]]

To write the above in the form w = N d we do

nodalsDegreOfFreedomVars={(Subscript[v, 1]),(Subscript[θ, 1]), (Subscript[v, 2]), (Subscript[θ, 2])};
{b,shapeFunctions}=Normal[CoefficientArrays[wNew,nodalsDegreOfFreedomVars]];
wNew=shapeFunctions.nodalsDegreOfFreedomVars;
Print["w(x)=",wNew]

Now that we have found the shape functions, we can find the B matrix. Notice that for beam, the strain matrix B is given by (((d^2)w)/(dx^2))

Β=D[shapeFunctions,{x,2}]

and now Π, the element potential energy is found (here I only consider strain energy, no loading is needed since we only need to find the element stiffness matrix, not solve a system for unknowns).

Ν={shapeFunctions};
Β={Β};
d=Transpose[{nodalsDegreOfFreedomVars}];
Dimensions[Ν]
Dimensions[d]
Dimensions[Β]

Π=(1/2)(∫_0^L)(Ε Ι Transpose[d].Transpose[Β].Β.d)\[DifferentialD]x;
Π=Π[[1,1]]

Now that we have the element potential energy, we can generate the element stiffness matrix

kmat=(∫_0^L)(Ε Ι Transpose[Β].Β)\[DifferentialD]x;
Print[ToString[((Ε Ι)/(L^3)),FormatType->TraditionalForm]<>ToString[kmat/(((Ε Ι)/(L^3))),FormatType->TraditionalForm]<>ToString[d,FormatType->TraditionalForm]]