Finite element analysis of a Variable length blade wind turbine: calculations of natural frequencies

%****************************************************************************
%****************************************************************************
%********************VARIBLADEANALYSIS **********************************
%****************************************************************************
%************************* Written by ******************************************
%****************************************************************************
%******************* TARTIBU KWANDA 2008**********************************
%****************************************************************************
%****************************************************************************
%****************************************************************************
%Variable length blade finite element program.  Two portions of beam: hollow beam %(fixed portion) and solid beam (moveable portion). Solves for eigenvalues and eigenvectors of %stepped beam with user defined dimensions and configurations. This program can calculate the %natural frequencies of different beam geometries and configurations ( position of moveable portion)
%default values are included in the program for the purpose of showing how to input data.
    echo off
    clf;
    clear all;
        inp = input('Input "1" to enter beam dimensions, "Enter" to use default ... ');
    if (isempty(inp))
        inp = 0;
    else
    end
   
    if inp == 0
%   input beam's geometries and material's properties
%   xl(i) = length of element (step)i
%   w(i) = width of element (step)i
%   t(i) = thickness of element (step)i
%   e = Young's modulus
%   bj = global degree of freedom number corresponding to the jth local degree
%   of freedom of element i
%   a(i) = area of cross section of element i
%   ne = number of elements
%   n = total number of degree of freedom
%   no = number of nodes
%   wa = width of hollow for the first portion of stepped beam
%   ta = thickness of hollow for the first portion of stepped beam
    format long
    beamconfiguration = 1;
    n1 = 10;
    n2 = 10;
    l1 = 1000;
    l2 = 1000;
    h1 = 1000/10;
    h2 = 1000/10;
    xl = [h1*ones(1,n1) h2*ones(1,beamconfiguration-1)];
    w1 = 60;
    w2 = 50;
    w = [w1*ones(1,n1) w2*ones(1,beamconfiguration-1)];
    t1 = 20;
    t2 = 10;
    t = [t1*ones(1,n1) t2*ones(1,beamconfiguration-1)];
    wh = 50;
    th = 10;
    wa = [wh*ones(1,beamconfiguration-1) 0*ones(1,n2)];
    ta = [th*ones(1,beamconfiguration-1) 0*ones(1,n2)];
%   configurations of stepped beam:
%   beamconfiguration ? n1 & n2
%   for configuration 1, w = [w1*ones(1,n1) w2*ones(1,1-1)],
%   t = [t1*ones(1,n1) t2*ones(1,1-1)]
%   wa = [wh*ones(1,1-1) 0*ones(1,n2)],
%   ta = [th*ones(1,1-1) 0*ones(1,n2)]
%   for configuration 2, w = [w1*ones(1,n1) w2*ones(1,2-1)]
%   t = [t1*ones(1,n1) t2*ones(1,2-1)]
%   wa = [wh*ones(1,2-1) 0*ones(1,n2)]
%   ta = [th*ones(1,2-1) 0*ones(1,n2)]
%   for configuration 3, w = [w1*ones(1,n1) w2*ones(1,3-1)]
%   t = [t1*ones(1,n1) t2*ones(1,3-1)]
%   wa = [wh*ones(1,3-1) 0*ones(1,n2)]
%   ta = [th*ones(1,3-1) 0*ones(1,n2)]
%   for configuration 4, w = [w1*ones(1,n1) w2*ones(1,4-1)]
%   t = [t1*ones(1,n1) t2*ones(1,4-1)]
%   wa = [wh*ones(1,4-1) 0*ones(1,n2)]
%   ta = [th*ones(1,4-1) 0*ones(1,n2)]
%   for configuration 5, w = [w1*ones(1,n1) w2*ones(1,5-1)]
%   t = [t1*ones(1,n1) t2*ones(1,5-1)]
%   wa = [wh*ones(1,5-1) 0*ones(1,n2)]
%   ta = [th*ones(1,5-1) 0*ones(1,n2)]
    xi = [40000 40000 40000 40000 40000 40000 40000 40000 40000 40000];
    a = [1200 1200 1200 1200 1200 1200 1200 1200 1200 1200];
    e = 206e+6;
    rho = 7.85e-6;
    bj = [1 2 3 4;3 4 5 6;5 6 7 8;7 8 9 10;9 10 11 12;11 12 13 14;13 14 15 16;15 16 17 18;17 18 19 20;19 20 21 22;21 22 23 24;23 24 25 26;25 26 27 28;27 28 29 30;29 30 31 32; 31 32 33 34;33 34 35 36;35 36 37 38;37 38 39 40;39 40 41 42];
    ne = 20;
    n = 22;
    else
       

    beamconfiguration = input ('Input beamconfiguration of stepped beam, default 1 ... ');
    if (isempty(beamconfiguration))
        beamconfiguration = 1;
    else
    end
   
    n1 = input ('Input number of elements for the first portion of stepped beam, default 10 ... ');
    if (isempty(n1))
        n1 = 10;
    else
    end
   
    n2 = input ('Input number of elements for the second portion of stepped beam, default 10 ... ');
    if (isempty(n2))
        n2 = 10;
    else
    end
   
    l1 = input ('Input length of first portion of stepped beam, default 1000, ... ');
    if (isempty(l1))
        l1 = 1000;
    else
    end
   
     l2 = input ('Input length of second portion of stepped beam, default 1000, ... ');
    if (isempty(l2))
        l2 = 1000;
    else
    end
   
    w1 = input ('Input width of first portion of stepped beam, default 60, ... ');
    if (isempty(w1))
            w1 = 60;
    else
    end
   
     w2 = input ('Input width of second portion of stepped beam, default 50, ... ');
    if (isempty(w2))
            w2 = 50;
    else
    end
   
    t1 = input ('Input thickness of first portion of stepped beam, default 20 , ... ');
    if (isempty(t1))
          t1 = 20;
    else
    end
   
    t2 = input ('Input thickness of second portion of stepped beam, default 10 , ... ');
    if (isempty(t2))
          t2 = 10;
    else
    end
   
    wh = input ('input width of hollow of first portion of stepped beam, default 50, ... ');
    if (isempty(wh))
        wh = 50;
    else
    end
   
    th = input ('input thickness of hollow of first portion of stepped beam, default 10, ... ');
    if (isempty(th))
        th = 10;
    else
    end
   
    e = input('Input modulus of material, mN/mm^2, default mild steel 206e+6 ... ');
    if (isempty(e))
        e=206e+6;
    else
    end
   
    rho = input('Input density of material,kg/mm^3 , default mild steel 7.85e-6 ... ');
    if (isempty(rho))
        rho = 7.85e-6;
    else
    end
   
    bj = input('Input global degree of freedom, default global degree of freedom [1 2 3 4;3 4 5 6;5 6 7 8;7 8 9 10;9 10 11 12;11 12 13 14;13 14 15 16;15 16 17 18;17 18 19 20;19 20 21 22;21 22 23 24;23 24 25 26;25 26 27 28;27 28 29 30;29 30 31 32; 31 32 33 34;33 34 35 36;35 36 37 38;37 38 39 40;39 40 41 42]; ... ');
    if (isempty(bj))
        bj =[1 2 3 4;3 4 5 6;5 6 7 8;7 8 9 10;9 10 11 12;11 12 13 14;13 14 15 16;15 16 17 18;17 18 19 20;19 20 21 22;21 22 23 24;23 24 25 26;25 26 27 28;27 28 29 30;29 30 31 32; 31 32 33 34;33 34 35 36;35 36 37 38;37 38 39 40;39 40 41 42];
    else
    end
    end
   

   
    ne = n1+beamconfiguration-1;
    h1 = l1/n1;
    h2 = l2/n2;
    xl = [h1*ones(1,n1) h2*ones(1,beamconfiguration-1)];
    w = [w1*ones(1,n1) w2*ones(1,beamconfiguration-1)];
    t = [t1*ones(1,n1) t2*ones(1,beamconfiguration-1)];
    wa = [wh*ones(1,beamconfiguration-1) 0*ones(1,n2)];
    ta = [th*ones(1,beamconfiguration-1) 0*ones(1,n2)];
   
%   calculate area (a), area moment of Inertia (xi) and mass per unit of length (xmas) of
%   the stepped beam
    a = (w.*t)-(wa.*ta);
%   define area moment of inertia according to flap-wise or edge-wise
%   vibration of the stepped beam.
    vibrationdirection = input('enter "1" for edge-wise vibration, "enter" for flap-wise vibration ... ');
    if (isempty(vibrationdirection))
        vibrationdirection = 0;
    else
    end
    if vibrationdirection == 0
        xi=((w.*t.^3)-(wa.*ta.^3))/12;
    else
        xi=((t.*w.^3)-(ta.*wa.^3))/12;
    end
   
    for i=1:ne
       xmas(i)=a(i)*rho;
    end
   
       

%   Size the stiffness and mass matrices
    no = ne+1;
    n = 2*no;
bigm = zeros(n,n);
bigk = zeros(n,n);
%   Now build up the global stiffness and consistent mass matrices, element
%   by element
for ii=1:ne
ai = zeros(4,n);
       i1=bj(ii,1);
       i2=bj(ii,2);
       i3=bj(ii,3);
       i4=bj(ii,4);
       ai(1,i1)=1;
       ai(2,i2)=1;
       ai(3,i3)=1;
       ai(4,i4)=1;
       xm(1,1)=156;
       xm(1,2)=22*xl(ii);
       xm(1,3)=54;
       xm(1,4)=-13*xl(ii);
       xm(2,2)=4*xl(ii)^2;
       xm(2,3)=13*xl(ii);
       xm(2,4)=-3*xl(ii)^2;
       xm(3,3)=156;
       xm(3,4)=-22*xl(ii);
       xm(4,4)=4*xl(ii)^2;
       xk(1,1)=12;
       xk(1,2)=6*xl(ii);
       xk(1,3)=-12;
       xk(1,4)=6*xl(ii);
       xk(2,2)=4*xl(ii)^2;
       xk(2,3)=-6*xl(ii);
       xk(2,4)=2*xl(ii)^2;
       xk(3,3)=12;
       xk(3,4)=-6*xl(ii);
       xk(4,4)=4*xl(ii)^2;
       for i=1:4
          for j=1:4
             xm(j,i)=xm(i,j);
             xk(j,i)=xk(i,j);
          end
       end
       for i=1:4
          for j=1:4
             xm(i,j)=(((xmas(ii)*xl(ii))/420))*xm(i,j);
             xk(i,j)=((e*xi(ii))/(xl(ii)^3))*xk(i,j);
          end
       end
       for i=1:n
          for j=1:4
             ait(i,j)=ai(j,i);
          end
       end
       xka=xk*ai;
       xma=xm*ai;
       aka=ait*xka;
       ama=ait*xma;
      for i=1:n
          for j=1:n
             bigm(i,j)=bigm(i,j)+ama(i,j);
             bigk(i,j)=bigk(i,j)+aka(i,j);
          end
      end
end
%   Application of boundary conditions
%   Rows and columns corresponding to zero displacements are deleted
   
     bigk(1:2,:) = [];
     bigk(:,1:2) = [];
       
     bigm(1:2,:) = [];
     bigm(:,1:2) = [];
   

%   Calculation of eigenvector and eigenvalue  
    [L, V] = eig (bigk,bigm)   
%   Natural frequency
    V1 = V.^(1/2)
    W = diag(V1)
     f=W/(2*pi)