Weird FEM issue (nat freqs)

  1. Hello,

    I recently used Matlab to find the natural frequencies of a clamped-clamped beam. It was fairly simple, as I construct the global mass and stiffness matrices. Then it's just a matter of using the eig() function. (For sake of simplicity I put beam properties all to 1.)

    When I choose how many elements (lets call it N) to use, it seems to give me natural frequencies which, seem to be smaller than usual. I somehow stumbled upon the fact that these frequencies are the literature values I have divided by N^2.

    Anyone got any ideas as why this is so? Here is my code. Thanks!! Edit: For example when N = 10, my first natural frequency is 0.2237 when the exact value is 22.37.

    Jeff

    Code (Text):
    syms y;
    syms z;
    area = 1;
    L = 1;
    p = 1;
    E = 1;
    Ig = 1;
    n = 10; %number of elements
    if n == 0
    stop
    end
    size = 4+((n-1)*2);
    %define global matrices
    Mg = zeros(size);
    Kg = zeros(size);
    %beam function
    A = [1 0 0 0 ; 0 1 0 0 ; 1 L L^2 L^3 ; 0 1 2*L 3*L^2];
    Ainv = inv(A);
    Aitrans = transpose(Ainv);
    Yh = [1 ; y ; y^2 ; y^3];
    Ya = [1 ; y];
    Ydotdot = diff(diff(Yh));


    %Solving for Mass matrix (kinetic energy)
    prod1 = Yh*transpose(Yh);
    M = p*area*Aitrans*int(prod1, y, 0, L)*Ainv;
    %solving for stiffness matrix (potential energy)
    K = (E*Ig)*Aitrans*int(Ydotdot*transpose(Ydotdot), y, 0, L)*Ainv;

    %creating the global mass matrix
    Mg(1:4, 1:4) = M(1:4, 1:4);
    Kg(1:4, 1:4) = K(1:4, 1:4);

    if n > 1
    i = 1;
    j = 1;
        for i=1:n-1
        Mg(j+2:j+5, j+2:j+5) = Mg(j+2:j+5, j+2:j+5) + M(1:4, 1:4);
        Kg(j+2:j+5, j+2:j+5) = Kg(j+2:j+5, j+2:j+5) + K(1:4, 1:4);
        j = j+2;
        end
       
    end

    %delete rows/columns based on clamped ends (boundary conditions)

    Mg(size,:) = [];
    Mg(size-1,:) = [];
    Mg(:,size) = [];
    Mg(:,size-1) = [];
    Mg(:,1) = [];
    Mg(:,1) = [];
    Mg(1,:) = [];
    Mg(1,:) = [];


    Kg(size,:) = [];
    Kg(size-1,:) = [];
    Kg(:,size) = [];
    Kg(:,size-1) = [];
    Kg(:,1) = [];
    Kg(:,1) = [];
    Kg(1,:) = [];
    Kg(1,:) = [];


    eigs = eig(Kg,Mg);
    sqrt(eigs)
     
  2. jcsd
  3. AlephZero

    AlephZero 7,298
    Science Advisor
    Homework Helper

    I haven't checked every line of the code, but I think you have each element of length 1, so the length of your beam depends on the number of elements in the model.

    If the length of the beam is [itex]l[/itex], the global mass is proportional to [itex]l[/itex] and the global stiffness proportional to [itex]1/l^3[/itex], so the frequences would be proportional to [itex]\sqrt{k/m} = 1/l^2[/itex].
     
  4. That was exactly it! Thank you kind sir!

    Jeff
     
Know someone interested in this topic? Share this thead via email, Google+, Twitter, or Facebook

Have something to add?