Parameters

a=1;
b=0.8;
E1=200 (10^9);
E2=150 (10^9);
ρ1=7850;
ρ2 =6000;
d1=0.04;
d2=0.03;
t1=0.05;
t2=0.06;
R1=0.1;
R2=0.08;

Visualization of the problem

Wall=Graphics[{Thick,Line[{{0,-(R1+R2)/2},{0,(R1+R2)/2}}]}];
Shaft1=Graphics[{EdgeForm[Thick],White,Rectangle[{0,-d1/2},{a,d1/2}]}];
Rotor1=Graphics[{EdgeForm[Thick],White,Rectangle[{a-t1/2,-R1},{a+t1/2,R1}]}];
Shaft2=Graphics[{EdgeForm[Thick],White,Rectangle[{a,-d2/2},{a+b,d2/2}]}];
Rotor2=Graphics[{EdgeForm[Thick],White,Rectangle[{a+b-t2/2,-R2},{a+b+t2/2,R2}]}];
Show[Wall,Shaft1,Shaft2,Rotor1,Rotor2]

Derived parameters

Moment of inertia of discs

m1=(R1^2) π t1 ρ1;
Θ1=m1 ((R1^2)/4+(t1^2)/12);
m2=(R2^2) π t2 ρ2;
Θ2=m2 ((R2^2)/4+(t2^2)/12);

Area moment of inertia of beam segments

I1=(d1^4) π/64;
I2=(d2^4) π/64;

Areas of cross sections

A1=(d1^2)π/4;
A2=(d2^2) π/4;

Computation of element matrices

Interpolation functions for bending (Hermitian polynomials)

v=Transpose[Table[{(Subscript[c, i])},{i,0,3}]].Table[{(ξ^i)},{i,0,3}];
φ=D[v,ξ];
Sv=Solve[{(v/.ξ->-1)[[1,1]]==v1,(φ/.ξ->-1)[[1,1]]==φ1,(v/.ξ->1)[[1,1]]==v2,(φ/.ξ->1)[[1,1]]==φ2},Table[(Subscript[c, i]),{i,0,3}]];
v=v[[1,1]]/.Sv;
Nd=Simplify[Coefficient[v[[1]],{v1,φ1,v2,φ2}]];
Nv=Table[If[Mod[i,2]==0,{(Le/2) Nd[[i]]},{Nd[[i]]}],{i,4}]

Plot of the Hermitian polynomials

Plot[Nd,{ξ,-1,1}]

Element stiffness matrix

Nv2=D[Nv,{ξ,2}];
Ke=Ie Ee (((2/Le))^3)(∫_(-1)^1)Nv2.Transpose[Nv2]\[DifferentialD]ξ;
Ke//MatrixForm

Element mass matrix

Me=ρe Ae Simplify[(Le/2)(∫_(-1)^1)Nv.Transpose[Nv]\[DifferentialD]ξ];
Me//MatrixForm

Calculation of global matrices

Particular element stiffness matrices

K1=Ke/.Le->a/.Ie->I1/.Ee->E1;
K2=Ke/.Le->b/.Ie->I2/.Ee->E2;

K1//MatrixForm
K2//MatrixForm

Particular element mass matrices

M1=Me/.Le->a/.Ae->A1/.ρe->ρ1;
M2=Me/.Le->b/.Ae->A2/.ρe->ρ2;

M1//MatrixForm
M2//MatrixForm

Construction of the global stiffness matrix

(Kg=ConstantArray[0,{6,6}])//MatrixForm

Kg[[1;;4,1;;4]]=Kg[[1;;4,1;;4]]+K1;
Kg//MatrixForm

Kg[[3;;6,3;;6]]=Kg[[3;;6,3;;6]]+K2;
Kg//MatrixForm

Construction of the global mass matrix of the beam

Mb=ConstantArray[0,{6,6}];
Mb[[1;;4,1;;4]]=Mb[[1;;4,1;;4]]+M1;
Mb[[3;;6,3;;6]]=Mb[[3;;6,3;;6]]+M2;
Mb//MatrixForm

Mass and inertia of rotors

Mr=ConstantArray[0,{6,6}];
Mr[[3,3]]=m1;
Mr[[4,4]]=Θ1;
Mr[[5,5]]=m2;
Mr[[6,6]]=Θ2;

Mr//MatrixForm

Global mass matrix

Mg=Mb+Mr;
Mg//MatrixForm

Calculation of eigenfrequencies

Constraints

All DoFs

CDoF={1,2};

All DoFs

ADoF = Table[i,{i,6}];
ADoF

Free DoFs

FDoF=Complement[ADoF,CDoF]

Condensed matrices

Condensed stiffness matrix

Kc=Kg[[FDoF,FDoF]];
Kc//MatrixForm

Condensed mass matrix

Mc=Mg[[FDoF,FDoF]];
Mc//MatrixForm

Solution of the frequency equation

Det[-α2 Mc+Kc]==0

Fr=Solve[Det[-α2 Mc+Kc]==0,α2]

(Sqrt[α2/.Fr])

It’s the same as the eigenvalue/vector problem of (M^(-1))K

Sort@Sqrt@Eigenvalues[Inverse[Mc].Kc]

Natural angular frequencies

Nfreq=(Sqrt[α2])/.Fr

In Hertz

Freq=Nfreq/(2π)

Mode shapes

Eigenvectors[Inverse[Mc].Kc]

Sort@Transpose@Eigensystem[Inverse[Mc].Kc]

A=(Sort@Transpose@Eigensystem[Inverse[Mc].Kc])[[All,2]]

U1=Table[{0,0,A[[i,1]],A[[i,2]]},{i,4}];
U2=Table[{A[[i,1]],A[[i,2]],A[[i,3]],A[[i,4]]},{i,4}];
v1d=U1.Nv/.Le->a;
v2d=U2.Nv/.Le->b;
v1=v1d/.ξ->(2 (x-0)/a-1);
v2=v2d/.ξ->(2 (x-a)/b-1);
Table[Show[Plot[v1[[j]],{x,0,a}],Plot[v2[[j]],{x,a,a+b}],PlotRange->All,ImageSize->{400,200}],{j,4}]

Manipulate[
(t^2),
{t,0,5,Appearance->"Labeled"}
]

Manipulate[
Show[Plot[v1[[j]]Cos[t],{x,0,a}],Plot[v2[[j]]Cos[t],{x,a,a+b}],PlotRange->{Automatic,0.6{-1,1}},ImageSize->{400,200}],
{{j,1,"Mode shape"},{1,2,3,4},Appearance->"Labeled"},
{t,0,2π,Appearance->"Labeled"}
]

Calculation of mode shapes

Aj=Table[Solve[((-α2/.Fr[[j]]) Mc+Kc).{{1},{Aj2},{Aj3},{Aj4}}=={{0},{0},{0},{0}},{Aj2,Aj3,Aj4}],{j,4}];
A=Table[{1,Aj2/.Aj[[j,1]],Aj3/.Aj[[j,1]],Aj4/.Aj[[j,1]]},{j,4}];

Element displacement vectors

U1=Table[{0,0,A[[i,1]],A[[i,2]]},{i,4}];
U2=Table[{A[[i,1]],A[[i,2]],A[[i,3]],A[[i,4]]},{i,4}];
v1d=U1.Nv/.Le->a;
v2d=U2.Nv/.Le->b;
v1=v1d/.ξ->(2 (x-0)/a-1);
v2=v2d/.ξ->(2 (x-a)/b-1);
Table[Show[Plot[v1[[j]],{x,0,a}],Plot[v2[[j]],{x,a,a+b}],PlotRange->All,ImageSize->{400,200}],{j,4}]