% FEAcomp. Solution of Example 4.1. 
%

%%% Compute the laminate constitutive matrices 
%%% and the equivalent Orthotropic laminate properties from
%%% LSS unidirectional laminae
%%%
%%% J.A. Mayugo, Febraury 2005
%%%

clear all;


n_file='FEAcomp_Ex401'
n_laminate='[(0/90)3]s'

%%% 3D properties, transversally isotropic material. 
%% Using material lamina axis (1 - long. fiber, 2 - transversal fiber)
E_1         = 145880;
E_2         = 13312;
G_12        = 4386;
G_23        = 4529;
nu_12       = 0.263;
nu_23       = 0.470;

%%% 3D Compliance and stiffness matrix in material directions
S_3d(1,1)   =1/E_1;
S_3d(2,2)   =1/E_2;
S_3d(3,3)   =S_3d(2,2);

S_3d(1,2)   =-nu_12/E_1;
S_3d(2,1)   =S_3d(1,2);
S_3d(1,3)   =S_3d(1,2);
S_3d(3,1)   =S_3d(1,3);

S_3d(2,3)   =-nu_23/E_2;
S_3d(3,2)   =S_3d(2,3);

S_3d(4,4)   =1/G_23;
S_3d(5,5)   =1/G_12;
S_3d(6,6)   =1/G_12;

Q_3d        =inv(S_3d)
 
%%% 2D Equivalent properties in 2D lamina

S_2d(1,1)   =1/E_1;
S_2d(2,2)   =1/E_2;
S_2d(1,2)   =-nu_12/E_1;
S_2d(2,1)   =S_2d(1,2);

S_2d(4,4)   =1/G_23;
S_2d(5,5)   =1/G_12;
S_2d(6,6)   =1/G_12;

Q_2d(1,1)=E_1/(1-nu_12^2*E_2/E_1);
Q_2d(1,2)=nu_12*E_2/(1-nu_12^2*E_2/E_1);
Q_2d(2,1)=Q_2d(1,2);
Q_2d(2,2)=E_2/(1-nu_12^2*E_2/E_1);
Q_2d(6,6)=G_12;
  
Q_2d(4,4)=G_23;
Q_2d(5,5)=G_12;

Q_2d

%%% The LSS 
% Laminate Input Data for transversaly isotropic materials, as
% [(0/90)2/+30/-30]s, angle (º), thickness (mm)

L(1,:)=[  0, 0.85];    % lamina #1
L(2,:)=[ 90, 0.85];    % lamina #2
L(3,:)=[  0, 0.85];    % lamina #3
L(4,:)=[ 90, 0.85];    % lamina #4
L(5,:)=[  0, 0.85];    % lamina #5
L(6,:)=[ 90, 0.85];    % lamina #6
L(7,:)=[ 90, 0.85];    % lamina #7
L(8,:)=[  0, 0.85];    % lamina #8
L(9,:)=[ 90, 0.85];    % lamina #9
L(10,:)=[  0, 0.85];    % lamina #10
L(11,:)=[ 90, 0.85];    % lamina #11
L(12,:)=[  0, 0.85];    % lamina #12

%% The first layer is the bottom layer....
%% Define equal number of L rows as laminae number
[M,N]=size(L);

%%% Compute Middle plane location
TH=0;
z=zeros(M,1);
th=zeros(M,1);
for i=1:M;
    TH=TH+L(i,2);
end;
th_mid=TH/2;

%%% Compute location and thickness of laminae
TH=0;
for i=1:M;
    th(i)=L(i,2);
    TH=TH+th(i);
    z(i)=TH-th(i)/2-th_mid;
end;

%%% Transform to laminate global axis (x - 0 degrees , y - 90 degrees)
T=zeros(6,6);
Ti=zeros(6,6);
Tg=zeros(6,6);
Q_eqv=zeros(6,6);
A=zeros(3,3);
B=zeros(3,3);
D=zeros(3,3);
H=zeros(2,2);
for i=1:M;
    theta = L(i,1)*pi/180;  % 
    m=cos(theta);
    n=sin(theta);
    a = [m  n 0;            % angle betwen x-y and 1-2
         -n  m 0;
         0  0 1];
     
% Transformation matrix T (for stresses) 
    for i1=1:3
        for j1=1:3
            if i1==j1 
                alpha=i1;
            else
                alpha=9-i1-j1;
            end
            for p1=1:3
                for q1=1:3
                    if p1==q1 
                        beta=p1;
                    else
                        beta=9-p1-q1;
                    end
                    if (alpha<=3)&(beta<=3)
                        T(alpha,beta)=a(i1,p1)*a(i1,p1);
                    elseif (alpha>3)&(beta<=3)
                        T(alpha,beta)=a(i1,p1)*a(j1,p1);
                    elseif (alpha<=3)&(beta>3)
                        T(alpha,beta)=a(i1,q1)*a(i1,p1)+a(i1,p1)*a(i1,q1);                        
                    elseif (alpha>3)&(beta>3)
                        T(alpha,beta)=a(i1,p1)*a(j1,q1)+a(i1,q1)*a(j1,p1);
                    end
                    if abs(T(alpha,beta))<1e-6 
                            T(alpha,beta)=0;
                    end
                end
            end
        end
    end
% Invert T
    Ti=inv(T);
% Transformation matrix Tg (for engineering strains)   
    Tg=Ti';    
% Stiffness an Compliance matrix in laminate global axis   
    Qb=Ti*Q_2d*Tg;
    Qb_3d=Ti*Q_3d*Tg;
%%% Compute 3d quasi-equivalent material  
   for j=1:6;
        for k=1:6;
% Compute A matrix
            Q_eqv(j,k)=Q_eqv(j,k)+(Qb_3d(j,k)*(th(i)))/TH;
        end;
    end;

%%% Compute A,B,D matrices 
   for j=1:3;
        for k=1:3;
            if k==3 
                k_=6;
            else
                k_=k;
            end;
            if j==3 
                j_=6;
            else
                j_=j;
            end;
% Compute A matrix
            A(j,k)=A(j,k)+Qb(j_,k_)*(th(i));
% Compute B matrix
            B(j,k)=B(j,k)+Qb(j_,k_)*(th(i)*z(i));
% Compute D matrix
            D(j,k)=D(j,k)+Qb(j_,k_)*(th(i)*z(i)^2+th(i)^3/12);
        end;
    end;
% Compute H matrix   
    for j=1:2;
        for k=1:2;
            j_=j+3;
            k_=k+3;
            H(j,k)=H(j,k)+(5/4)*Qb(j_,k_)*(th(i)-4/TH^2*(th(i)*z(i)^2+th(i)^3/12));
        end;
    end;    
    
end;

S_eqv=inv(Q_eqv);

E_xeqv=1/S_eqv(1,1);
E_yeqv=1/S_eqv(2,2);
E_zeqv=1/S_eqv(3,3);

nu_xyeqv=-E_xeqv*S_eqv(1,2);
nu_xzeqv=-E_xeqv*S_eqv(1,3);
nu_yzeqv=-E_yeqv*S_eqv(2,3);

G_xyeqv=1/S_eqv(6,6);
G_xzeqv=1/S_eqv(5,5);
G_yzeqv=1/S_eqv(4,4);

%%% OUTPUT DATA
format short e
title0=['Laminate' n_laminate]
titleA=['Casi equivalent stiffness 3D']

Q_eqv

title1=['A, B, D and H matrices']

ABD =[A B;B D]

H 

titleB=['Input data for ANSYS SHELL93 equivalent laminate'];

Exeqv = num2str(E_xeqv, '%5.0f');
Eyeqv = num2str(E_yeqv, '%5.0f');
Ezeqv = num2str(E_zeqv, '%5.0f');

Gxyeqv = num2str(G_xyeqv, '%5.0f');
Gxzeqv = num2str(G_xzeqv, '%5.0f');
Gyzeqv = num2str(G_yzeqv, '%5.0f');

nuxyeqv = num2str(nu_xyeqv, '%3.3f');
nuxzeqv = num2str(nu_xzeqv, '%3.3f');
nuyzeqv = num2str(nu_yzeqv, '%3.3f');

tline1a = ['uimp,1,ex,ey,ez,' Exeqv ',' Eyeqv ',' Ezeqv];
tline1b = ['uimp,1,gxy,gyz,gxz,' Gxyeqv ',' Gyzeqv ',' Gxzeqv];
tline1c = ['uimp,1,prxy,pryz,prxz,' nuxyeqv ',' nuyzeqv ',' nuxzeqv];

title2=['Input data for ANSYS SHELL99 with KEYOPTION(2)=2'];

A11 = num2str(A(1,1), '%6g');
A12 = num2str(A(1,2), '%6g');
A16 = num2str(A(1,3), '%6g');
A22 = num2str(A(2,2), '%6g');
A26 = num2str(A(2,3), '%6g');
A66 = num2str(A(3,3), '%6g');
H44 = num2str(H(1,1), '%6g');
H45 = num2str(H(1,2), '%6g');
H55 = num2str(H(2,2), '%6g');
B11 = num2str(B(1,1), '%6g');
B12 = num2str(B(1,2), '%6g');
B16 = num2str(B(1,3), '%6g');
B22 = num2str(B(2,2), '%6g');
B26 = num2str(B(2,3), '%6g');
B66 = num2str(B(3,3), '%6g');
D11 = num2str(D(1,1), '%6g');
D12 = num2str(D(1,2), '%6g');
D16 = num2str(D(1,3), '%6g');
D22 = num2str(D(2,2), '%6g');
D26 = num2str(D(2,3), '%6g');
D66 = num2str(D(3,3), '%6g');
thick = num2str(TH, '%6g');

tline1 = ['R,1,' A11 ',' A12 ',0,' A16 ',0,0'           '         ! real set #1, A11,A12,0,A16,0,0'];
tline2 = ['RMODIF,1,7,' A22 ',0,' A26 ',0,0'            '         ! real set #1, A22,0,A26,0,0'];
tline3 = ['RMODIF,1,16,' A66 ',0,0'                     '         ! real set #1, A66,0,0'];
tline4 = ['RMODIF,1,19,' H44 ',' H45                    '         ! real set #1, H44,H45'];
tline5 = ['RMODIF,1,21,' H55                            '         ! real set #1, H55'];
tline6 = ['RMODIF,1,22,' B11 ',' B12 ',0,' B16 ',0,0'   '         ! real set #1, B11,B12,0,B16,0,0'];
tline7 = ['RMODIF,1,28,' B22 ',0,' B26 ',0,0'           '         ! real set #1, B22,0,B26,0,0'];
tline8 = ['RMODIF,1,37,' B66 ',0,0'                     '         ! real set #1, B66,0,0'];
tline9 = ['RMODIF,1,43,' D11 ',' D12 ',0,' D16 ',0,0'   '         ! real set #1, D11,D12,0,D16,0,0'];
tline10 =['RMODIF,1,49,' D22 ',0,' D26 ',0,0'           '         ! real set #1, D22,0,D26,0,0'];
tline11 =['RMODIF,1,58,' D66 ',0,0'                     '         ! real set #1, D66,0,0'];
tline12 =['RMODIF,1,76,,' thick '          ! real set #1, Average density and thick'];

% disp(tline1);
% disp(tline2);
% disp(tline3);
% disp(tline4);
% disp(tline5);
% disp(tline6);
% disp(tline7);
% disp(tline8);
% disp(tline9);
% disp(tline10);
% disp(tline11);
% disp(tline12);

fid0 = fopen([n_file,'_ABDH.txt'],'w');

fprintf(fid0,'%10s\n\n',title0);

fprintf(fid0,'%10s\n\n',title1);
fprintf(fid0,'%12.4e\t %12.4e\t %12.4e\t %12.4e\t %12.4e\t %12.4e\n',ABD);

fprintf(fid0,'%10s\n\n',' ');
fprintf(fid0,'%12.4e\t %12.4e\n',H);
fprintf(fid0,'%10s\n\n',' ');

fprintf(fid0,'\n%10s\n\n',titleB);
fprintf(fid0,'%10s\n',tline1a);
fprintf(fid0,'%10s\n',tline1b);
fprintf(fid0,'%10s\n',tline1c);

fprintf(fid0,'\n%10s\n\n',title2);
fprintf(fid0,'%10s\n',tline1);
fprintf(fid0,'%10s\n',tline2);
fprintf(fid0,'%10s\n',tline3);
fprintf(fid0,'%10s\n',tline4);
fprintf(fid0,'%10s\n',tline5);
fprintf(fid0,'%10s\n',tline6);
fprintf(fid0,'%10s\n',tline7);
fprintf(fid0,'%10s\n',tline8);
fprintf(fid0,'%10s\n',tline9);
fprintf(fid0,'%10s\n',tline10);
fprintf(fid0,'%10s\n',tline11);
fprintf(fid0,'%10s\n',tline12);

fclose(fid0);
