Code covered by the BSD License  

Highlights from
Co-Blade: Software for Analysis and Design of Composite Blades

image thumbnail

Co-Blade: Software for Analysis and Design of Composite Blades

by

 

18 Sep 2012 (Updated )

Analysis and design of composite blades for wind and hydrokinetic turbines

cltAnalysis(lamNum,nPlies,tPly,fibAng,E11,E22,G12,nu12,density)
function CLT = cltAnalysis(lamNum,nPlies,tPly,fibAng,E11,E22,G12,nu12,density)
    
d2r = pi/180;   % convert degrees to radians

%% Compute total thickness of the laminate
t = sum( nPlies.*tPly );

%% Compute mass per unit volume
massByV = sum(nPlies.*tPly.*density) / t; % (kg/m^3)

%% Compute ply (lamina) interface locations
N    = length(lamNum);
z    = zeros(N+1,1);
z(1) = -t/2;
for k = 1:N
    z(k+1) = z(k) + nPlies(k)*tPly(k);
end

%% Calculate the reduced stiffness matrix, Q, for each lamina
Q = zeros(3,3,N);
for k = 1:N    
    Q(1,1,k) = E11(k)^2              / (E11(k) - nu12(k)^2 * E22(k));
    Q(2,2,k) = E11(k)*E22(k)         / (E11(k) - nu12(k)^2 * E22(k));
    Q(1,2,k) = nu12(k)*E11(k)*E22(k) / (E11(k) - nu12(k)^2 * E22(k));
    Q(2,1,k) = Q(1,2,k);
    Q(3,3,k) = G12(k);     
end 

%% Calculate the transformed reduced stiffness matrix, Qbar, and transformation matrix, T, for each lamina
Qbar = zeros(3,3,N);
T    = zeros(3,3,N);
for k = 1:N
    c = cos( fibAng(k)*d2r );
    s = sin( fibAng(k)*d2r );
    
    Qbar(1,1,k) = (c^4)*Q(1,1,k) + 2*(c^2)*(s^2)*(Q(1,2,k) + 2*Q(3,3,k)) + (s^4)*Q(2,2,k);
    Qbar(1,2,k) = (c^4 + s^4)*Q(1,2,k) + (c^2)*(s^2)*(Q(1,1,k) + Q(2,2,k) - 4*Q(3,3,k));
    Qbar(2,1,k) = Qbar(1,2,k);
    Qbar(1,3,k) = (c^3)*s*(Q(1,1,k) - Q(1,2,k) - 2*Q(3,3,k)) - c*(s^3)*(Q(2,2,k) - Q(1,2,k) - 2*Q(3,3,k));
    Qbar(3,1,k) = Qbar(1,3,k);
    Qbar(2,2,k) = (s^4)*Q(1,1,k) + 2*(c^2)*(s^2)*(Q(1,2,k) + 2*Q(3,3,k)) + (c^4)*Q(2,2,k);
    Qbar(2,3,k) = c*(s^3)*(Q(1,1,k) - Q(1,2,k) - 2*Q(3,3,k)) - s*(c^3)*(Q(2,2,k) - Q(1,2,k) - 2*Q(3,3,k));
    Qbar(3,2,k) = Qbar(2,3,k);
    Qbar(3,3,k) = (c^2)*(s^2)*(Q(1,1,k) + Q(2,2,k) - 2*Q(1,2,k) - 2*Q(3,3,k)) + (c^4 + s^4)*Q(3,3,k);
    
    T(1,1,k) = c^2;
    T(1,2,k) = s^2;
    T(1,3,k) = 2*c*s;
    T(2,1,k) = s^2;
    T(2,2,k) = c^2;
    T(2,3,k) = -2*c*s;
    T(3,1,k) = -c*s;
    T(3,2,k) = c*s;
    T(3,3,k) = c^2 - s^2;
 
end 

%% Calculate the A, B, D, and ABD matrices
A = zeros(3,3); % extensional stiffness matrix
B = zeros(3,3); % coupling matrix
D = zeros(3,3); % bending stiffness matrix
for i = 1:3
    for j = 1:3
        sumA = 0;
        sumB = 0;
        sumD = 0;
        for k = 1:N
            sumA = sumA + Qbar(i,j,k)*(z(k+1) - z(k));
            sumB = sumB + Qbar(i,j,k)*(z(k+1)^2 - z(k)^2);
            sumD = sumD + Qbar(i,j,k)*(z(k+1)^3 - z(k)^3);
        end
        A(i,j) = sumA;
        B(i,j) = (1/2) * sumB;
        D(i,j) = (1/3) * sumD;     
    end
end
ABD = [A B; 
       B D];

%% The effective elastic properties of the laminate
a      = inv(A);
E_eff  = 1/(t*a(1,1));   % effective Young's modulus of the laminate
G_eff  = 1/(t*a(3,3));   % effective shear modulus of the laminate
nu_eff = -a(1,2)/a(1,1); % effective Poisson ratio of the laminate

%% Collect output
CLT.t       = t;
CLT.E_eff   = E_eff;
CLT.G_eff   = G_eff;
CLT.nu_eff  = nu_eff;
CLT.massByV = massByV;
CLT.A       = A;
CLT.B       = B;
CLT.D       = D;
CLT.ABD     = ABD;
CLT.Qbar    = Qbar;
CLT.T       = T;
CLT.z_i     = z;

end %function cltAnalysis

Contact us