Plotting eigenmodes in 3D

Hi!

I have a code that calculates eigenfrequencies for a frame structure thruogh the "eig"-command. I have the stiffness and massmatrix as input and get the eigenvalues and eigenvectors as output. What I want to do is to plot a 3D-plot of the eigenmodes of the structure. I want to see which frequency that corresponds to bending, torsion ...

I have not found any good ways to plot the eigenvector in 3D, any advise?

Thank you!

[Van Tien Nguyen]

Hi, I am looking for a code to solve eigenvalue problem for 3D frame building. Can you exchange your code? Thank you in advance.

Answer 1

[KSSV]

This shoul dhelp: https://in.mathworks.com/matlabcentral/fileexchange/30970-natural-frequencies-buckling-loads-of-columns-using-finite-element-method

[philip singh]

Thanks for the quick answer.

The beam in this code is unfortunately a 1D-beam. I have a 3D-frame, where the stiffnesselementmatrix is 12x12 matrix and the assembled is 1356x1356. Therfore the plot-function in the link is not applicable on my frame.

Answer 2

[Christine Tobler]

Depending on the format of your 3D-frame, maybe you could use pdeplot3d.

[philip singh]

% Code to generate inputfile
close all
clear all, close all, clc
format long
L = 5.6;
le = 0.05;
x=0:le:L;
x=x';
Nn = 226;
Ne = 228;
nodeCoordinates=zeros(length(x),3);

% node-coordinates
for i=1:length(x)
    nodeCoordinates(i,1)=x(i);
    nodeCoordinates(i,2)=0;
    nodeCoordinates(i,3)=0;
end
x=x';
x=fliplr(x);
for i=1:length(x)
    nodeCoordinates(length(x)+i,1)=x(i);
    nodeCoordinates(length(x)+i,2)=0.76;
    nodeCoordinates(length(x)+i,3)=0;
end
A1 = repelem(1:1:Nn,1); 
node = [A1;nodeCoordinates(:,1)';nodeCoordinates(:,2)';nodeCoordinates(:,3)'];



vektor1=zeros(Ne,1);
pos_TB_V1=1; % pos. cross-member left 1
pos_TB_V2=26; % pos. cross-member left 2
pos_TB_V3=52; % pos. cross-member left 3
pos_TB_V4=70; % pos. cross-member left 4
pos_TB_V5=110; % pos. cross-member left 5

% element position
for i=1:Ne
    
    if i == pos_TB_V1
        vektor1(i,1)=pos_TB_V1;
        vektor1(i+1,1)=pos_TB_V1;
    
    elseif i==pos_TB_V2
        vektor1(i+1,1)=pos_TB_V2;
        vektor1(i+2,1)=pos_TB_V2;
     
    elseif i == pos_TB_V3
        vektor1(i+2,1)=pos_TB_V3;
        vektor1(i+3,1)=pos_TB_V3;
     
    elseif i == pos_TB_V4
        vektor1(i+3,1)=pos_TB_V4;
        vektor1(i+4,1)=pos_TB_V4;
        
    elseif i == pos_TB_V5
        vektor1(i+4,1)=pos_TB_V5;
        vektor1(i+5,1)=pos_TB_V5;
    end
    
end

vektor1(3:pos_TB_V2,1)=2:pos_TB_V2-1;
vektor1(pos_TB_V2+3:pos_TB_V3+1,1)=pos_TB_V2+1:pos_TB_V3-1;
vektor1(pos_TB_V3+4:pos_TB_V4+2,1)=pos_TB_V3+1:pos_TB_V4-1;
vektor1(pos_TB_V4+5:pos_TB_V5+3,1)=pos_TB_V4+1:pos_TB_V5-1;
vektor1(pos_TB_V5+6:Ne,1)=pos_TB_V5+1:Ne-5;
vektor1((Nn+10)/2)=[];
vektor1(Ne,1)=Nn-2;

vektor2=zeros(Ne,1); 
pos_TB_H1=Nn-pos_TB_V1+1; % pos. cross-member right 1
pos_TB_H2=Nn-pos_TB_V2+1; % pos. cross-member right 2
pos_TB_H3=Nn-pos_TB_V3+1; % pos. cross-member right 3
pos_TB_H4=Nn-pos_TB_V4+1; % pos. cross-member right 4
pos_TB_H5=Nn-pos_TB_V5+1; % pos. cross-member right 5


for i=1:Ne+1
    
    if i == pos_TB_V1
        vektor2(i,1)=pos_TB_H1;
        
    elseif i==pos_TB_V2
        vektor2(i,1)=pos_TB_V2;
        vektor2(i+1,1)=pos_TB_H2;
        
    elseif i == pos_TB_V3
        vektor2(i+1,1)=pos_TB_V3;
        vektor2(i+2,1)=pos_TB_H3;
        
    elseif i == pos_TB_V4
        vektor2(i+2,1)=pos_TB_V4;
        vektor2(i+3,1)=pos_TB_H4;
        
    elseif i == pos_TB_V5
        vektor2(i+3,1)=pos_TB_V5;
        vektor2(i+4,1)=pos_TB_H5;
    end
end

vektor2(2:pos_TB_V2-1,1)=2:pos_TB_V2-1;
vektor2(pos_TB_V2+2:pos_TB_V3,1)=pos_TB_V2+1:pos_TB_V3-1;
vektor2(pos_TB_V3+3:pos_TB_V4+1,1)=pos_TB_V3+1:pos_TB_V4-1;
vektor2(pos_TB_V4+4:pos_TB_V5+2,1)=pos_TB_V4+1:pos_TB_V5-1;
vektor2(pos_TB_V5+5:Ne,1)=pos_TB_V5+1:Ne-4;
vektor2((Ne+8)/2)=[];
vektor2(Ne,1)=Ne-3;
length(vektor2);
vektor2;

E = 2.1E11;
A = 3.52E-3;
Iz = 2.414E-6;
Iy = 3.758E-6;
G = 8.077E10;
K = 7.51E-8;
Rho = 7.85E3;

% element parameters
for i = 1:Ne
   
    vektorE(i,1)= E;
    vektorA(i,1)= A;
    vektorIz(i,1)= Iz;
    vektorIy(i,1)= Iy;
    vektorG(i,1)= G;
    vektorK(i,1)= K;
    vektorRho(i,1)= Rho;
end

B1 = repelem(1:1:Ne,1); 
element = [B1;vektor1';vektor2';vektorE';vektorA';vektorIz';vektorIy';vektorG';vektorK';vektorRho'];

% wrights all data to inputfile
fid = fopen('Ram3D.inp','w');
fprintf(fid,'%s\n','*NODE [Node #, X, Y, Z]');
fprintf(fid,'%3.0f, %6.2f, %6.2f, %6.2f,\n',node);

fprintf(fid,'%s\n','*ELEMENT [El #, 1st Node #, 2nd Node #, E, A, Iz, Iy, G, J, rho]');
fprintf(fid,'%3.0f, %6.0f, %6.0f, %15.0f, %15.10f, %15.10f, %15.10f, %15.0f, %15.10f, %6.0f,\n',element);
fprintf(fid,'%s','*END');
fclose(fid);

%type Ram3D.inp

% Code for calculating eigenfrequencies

% Variables:
%    Problem size
%       Nn:         total number of nodes
%       Ne:         total number of elements
%       No_of_Eqns: total number of equations (degrees of freedom)
%    Vectors & Matrices:
%      node:    Columns: 1 = node #, 2 = x-coord
%      element: Columns: 1 = 1st node #, 2 = 2nd node #, 3 = E modulus,
%                        4 = Area, 5 = moment of intertia Iz, 
%                        6 = moment of intertia Iy, 7 = G, 8 = K, 9 = Rho
%------------------------------------------------------------
% 1. Read in-put data from file = frame2D.inp
%       The global force vector is established here!
fid = fopen('Ram3D.inp','r');
[node,element,Nn,Ne] = read_input(fid);
fclose(fid);
%------------------------------------------------------------
% 2. Assembly of (i) global stiffness matrix of the frame
%       and (ii) global force vector.
%       Loop over all elements (compute the element stiffness matrix
%       for each element in the global coordinate system and add it
%       to the system stiffness matrix).
%       Define the dimension of the matrix and set it to zero.
%
%       Element stiffnes matrix in global coordinate system
%           Ke(i,j)  dimension: 6x6 
%       System stiffnes matrix in global coordinate system
%           Ksystem(i,j)  dimension: No_of_Eqns x No_of_Eqns
%
No_of_Eqns = 6*Nn;
Ksystem = zeros(No_of_Eqns,No_of_Eqns);
Msystem = zeros(No_of_Eqns,No_of_Eqns);

for e = 1 : Ne
    Enode1 = round(element(e,1)); Enode2 = round(element(e,2));
    Eemod = element(e,3); Earea = element(e,4); EIz = element(e,5);...
    EIy = element(e,6); EG = element(e,7); EK = element(e,8); Erho = element(e,9);
    [Ke,Me] = Matrix(Eemod,Earea,EIz,EIy,EG,EK,Erho,node(Enode1,:),node(Enode2,:));
    ekv(6) = 6*Enode1;
    ekv(5) = ekv(6) - 1;
    ekv(4) = ekv(5) - 1;
    ekv(3) = ekv(4) - 1;
    ekv(2) = ekv(3) - 1;
    ekv(1) = ekv(2) - 1;
    ekv(12) = 6*Enode2;
    ekv(11) = ekv(12) - 1;
    ekv(10) = ekv(11) - 1;
    ekv(9) = ekv(10) - 1;
    ekv(8) = ekv(9) - 1;
    ekv(7) = ekv(8) - 1;
    for i = 1 : 12
        k = ekv(i);
        for j = 1 : 12
            l = ekv(j);
            Ksystem(k,l) = Ksystem(k,l) + Ke(i,j);
            Msystem(k,l) = Msystem(k,l) + Me(i,j);
        end
    end
end

% The matrix is badly conditioned,
% adding a small element to the diagonal will help,
c=7.96e-6;
K_G = Ksystem + c*eye(size(Ksystem));
Ksystem=K_G;

[V,W]=eig(Ksystem,Msystem); 

omega = sqrt(diag(W));
f=(omega/(2*pi));
F_mat = f(1:10);

disp('Egenfrekvenserna för systemet:')

for i=1:length(F_mat)
    disp(['f',num2str(i),' = ',num2str(F_mat(i)),' Hz'])
end

%*************************************************************
%           FUNCTIONS     and     SUBROUTINES
%*************************************************************

% --------Compute the elementstiffness & mass matrix ----------

function [Ke, Me] = Matrix(E,A,Iz,Iy,G,K,rho,node1,node2)

%  Retrieve coordinates of element nodes
x1 = node1(1);  x2 = node2(1);  y1 = node1(2);  y2 = node2(2);
z1 = node1(3);  z2 = node2(3);
le = sqrt((x2-x1)^2 + (y2-y1)^2 + (z2-z1)^2);  l1 = le;

Ke=zeros(12,12);
Ke(1,1)  = E*A/l1;
Ke(1,7)  = -Ke(1,1);
Ke(2,2)  = 12*E*Iz/(l1*l1*l1);
Ke(2,6)  = 6*E*Iz/(l1*l1);
Ke(2,8)  = -Ke(2,2);
Ke(2,12) = Ke(2,6);
Ke(3,3)  = 12*E*Iy/(l1*l1*l1);
Ke(3,5)  = -6*E*Iy/(l1*l1);
Ke(3,9)  = -Ke(3,3);
Ke(3,11) = Ke(3,5);
Ke(4,4)  = G*K/l1;
Ke(4,10) = -Ke(4,4);
Ke(5,3) = Ke(3,5);
Ke(5,5) = 4*E*Iy/l1;
Ke(5,9) = -Ke(3,5);
Ke(5,11) = 2*E*Iy/l1;
Ke(6,2) = Ke(2,6);
Ke(6,6) = 4*E*Iz/l1;
Ke(6,8) = -Ke(2,6);
Ke(6,12) = 2*E*Iz/l1;
Ke(7,1)  = -Ke(1,1);
Ke(7,7) = Ke(1,1);
Ke(8,2) = -Ke(2,2);
Ke(8,6) = -Ke(2,6);
Ke(8,8) = Ke(2,2);
Ke(8,12) = -Ke(2,6);
Ke(9,3) = -Ke(3,3);
Ke(9,5) = -Ke(3,5);
Ke(9,9) = Ke(3,3);
Ke(9,11) = -Ke(3,5);
Ke(10,4) = -Ke(4,4);
Ke(10,10) = Ke(4,4);
Ke(11,3) = Ke(3,5);
Ke(11,5) = Ke(5,11);
Ke(11,9) = -Ke(3,5);
Ke(11,11) = Ke(5,5);
Ke(12,2) = Ke(2,6);
Ke(12,6) = Ke(6,12);
Ke(12,8) = -Ke(2,6);
Ke(12,12) = Ke(6,6);


Me=zeros(12,12);
Me(1,1) = rho*A*l1/2;
Me(2,2) = Me(1,1);
Me(3,3) = Me(1,1);
Me(4,4) = rho*(Iy+Iz)*l1/2;
Me(5,5) = rho*A*l1.^3/24;
Me(6,6) = Me(5,5);
Me(7,7) = Me(1,1);
Me(8,8) = Me(1,1);
Me(9,9) = Me(1,1);
Me(10,10) = Me(4,4);
Me(11,11) = Me(5,5);
Me(12,12) = Me(5,5);
% ------ Coordinate transformation (From local to global system) -------
if x1 == x2 && y1 == y2
if z2>z1

    Lambda = [0 0 1 ;0 1 0 ;-1 0 0];
else
Lambda = [0 0 -1 ;0 1 0 ;1 0 0];
end
else
CXx = (x2-x1)/l1;
CYx = (y2-y1)/l1;
CZx = (z2-z1)/l1;
D = sqrt(CXx*CXx + CYx*CYx);
CXy = -CYx/D;
CYy = CXx/D;
CZy = -CZx/D;
CXz = (-CXx*CZx)/D;
CYz = (-CYx*CZx)/D;
CZz = D;

Lambda = [CXx CYx CZx ;CXy CYy CZy ;CXz CYz CZz];
end

R = [Lambda zeros(3,9); zeros(3) Lambda zeros(3,6);
zeros(3,6) Lambda zeros(3);zeros(3,9) Lambda];

Ke = R'*Ke*R;
Me = R'*Me*R;

end
%*************************************************************
% This routine reads in-put data from the file "Ram3D.inp"
function [node,element,Nn,Ne] = read_input(fid)
data = [];
% Node numbers and coordinates
data = str2num(fgetl(fid));
data = str2num(fgetl(fid));
while (length(data) > 0)
    Nnum = round(data(1));
    node(Nnum,1:3) = data(2:4);
    data = str2num(fgetl(fid));
end
Nn = length(node(:,1));

% Elements (connectivity)
data = str2num(fgetl(fid));
while (length(data) > 0)
    Enum = round(data(1));
    element( Enum, 1:9) = data(2:10);
    data = str2num(fgetl(fid));
end
Ne = length(element(:,1));

end

This is my code. The first part creats a inputfile that is later used. The inputfile contains the nodecoordinates and the elementparameters. If you run the code the first ten frequencies is printed in the command window, it is for these ten frequencies that I want to plot corresponding eigenmodes (as a 3D-figure). I want to see how the eigenmodes looks like. I have tried to plot the eigenvectors with different commands but have not succeeded in getting a good plot.

[Christine Tobler]

Sorry I can't quite understand from your code what kind of elements you are modeling. Can you describe what geometry this is based on. Where are the coordinates stored? Is this a tetrahedral mesh, or rectangular box mesh, or something else? Where are these elements defined?

[philip singh]

It's my fault, I should have commented the code more but thanks for the help.

The frame consists out of two long beams parallel to each other and five shorter cross-members between the two long beams. Each element is 5 centimeter and the length of the total frame is 5.6 meters. The elements and coordinates are defined in the inputfile that is created in the beginning of the code above. The coordinates are stored in a variable called "node" and the elements with its area, stiffness, density... is stored in a variable called "element".

Does this answer your question?

[Christine Tobler]

I'm still a bit confused: the third column of node seems to be all zero, so is this really a 2D problem? That would be much simpler to plot.

Also, the element matrix seems to contain only two nodes, so is each element simply a line? I would have expected a triangle / rectangle / tetrahedron / ... for a 3D problem.

Anyway, here are some functions that may be helpful to visualize PDe solutions in MATLAB: trisurf, tetramesh, pdeplot, pdeplot3d (particularly the nodes, elements syntax).

Answer 3

[jice zeng]

Hi，do you figure out your problem? I also have a 3D frame structure to be analyzed for its dynamic properties, such as frequency and mode shape, I have trouble in plotting 3d mode shape, so can you share your idea to plot mode shape with me?

[Ravi Kumar]

Check this example in PDE Toolbox:

https://www.mathworks.com/help/pde/examples/vibration-of-a-square-plate.html