% Elastic micromechanics, transversely isotropic fibers [Barbero & Luciano 1995].
%  Copyright: Ever J. Barbero
%             Department of Mechanical and Aerospace Engineering
%             West Virginia University, Morgantown, WV 26505-6106
%  28 Feb 2006, MATLAB 7 
%  IF YOU PUBLISH WORK BENEFITING FROM THIS M-FILE, PLEASE CITE IT AS:
%  E. J. Barbero and R. Luciano, "Micromechanical Formulas for the Relaxation Tensor of 
%  Linear Viscoelastic Composites with Transversely Isotropic Fibers," Int. J. of Solids 
%  Structures, 32(13), 1859-1872, (1995).
%  http://barbero.cadec-online.com/feacm/Examples/Ch6ex/PMMTE.m
clear all; close all;
format short g
time = [.01;.1;1;2;3;5;10;20;30;50;100;200;300];
syms s complex;
disp('Periodic Microstructure with transversely isotropic fibers and elastic matrix');
disp('=============================================================================');
%% Matrix Vinylester 
% elastic solution
E_0=0.592; nu_m = 0.311; E_m = E_0;

%Fiber: transversely isotropic																					
E_a	= 24.42; E_t = 3.6; nu_a = 0.443; nu_t = 0.005; G_a = 6.4;
G_t = E_t/2/(1+nu_t); % calculated
V_f = 0.6;	

% calculation section
C16=E_m;
D16=nu_m;
F15=V_f;
B12=E_a;
C12=E_t;
D12=nu_a;
E12=nu_t;
F12=G_a;
G12=G_t;

% restrictions
H12=1-(D12*D12*C12/B12);	
I12=1-(D12*D12*C12/B12);	
J12=1-(E12*E12);	
K12=1-(D12*D12*C12/B12)-(E12*E12)-(D12*D12*C12/B12)-(2*(D12*C12/B12)*(E12)*(D12));
if ((H12<=0) | (I12<=0) | (J12<=0) | (K12<=0))
    disp('Error: restriction on FIBER elastic constants')
end

% matrix
E35=C16*D16/((1+D16)*(1-2*D16));
F35=C16/(2*(1+D16));	
G35=D16;
% fiber
E38=J12*B12/K12;	
F38=(D12+E12*D12)*C12/K12;	
G38=I12*C12/K12;
H38=(E12+((C12*D12/B12)*D12))*C12/K12;	
I38=F12;
% S terms
D42=0.49247-0.47603*F15-0.02748*F15^2;	
E42=0.36844-0.14944*F15-0.27152*F15^2;	
F42=0.12346-0.32035*F15+0.23517*F15^2;
% a1 a2 a3 a4 etc
J34=(4*F35^2)-(2*F35*G38)+(6*E35*F35)-(2*E38*F35)-(2*F35*H38)+(H38*E38)+(4*E35*F38)-(2*F38^2)-(E35*G38)-(2*E38*E35)+(E38*G38)-(E35*H38);	
K34=(8*F35^3)-(8*F35*F35*G38)+(12*F35*F35*E35)-(4*F35*F35*E38)-(2*F35*H38*H38)+(4*F35*E35*H38)+(4*F35*E38*G38)-(8*F35*E35*G38)-(4*F35*F38*F38)+(2*F35*G38*G38)-(4*F35*E38*E35)+(8*F35*E35*F38)+(2*F35*E38*G38)+(4*F38*H38*E35)-(4*F38*G38*E35)-(2*E35*E38*H38)-(2*H38*F38*F38)+(H38*H38*E38)+(2*G38*F38*F38)-(E38*G38*G38)+(E35*G38*G38)-(E35*H38*H38);	
L34=(((4*F35^2)+(4*F35*E35)-(2*E38*F35)-(2*F35*G38)-(E38*E35)-(E35*G38)-(F38^2)+(E38*G38)+(2*E35*F38))/(K34))-(((D42)-(E42/(2-2*G35)))/(F35));	
M34=(-((-2*F35*H38)+(2*E35*F35)-(E35*H38)-(E38*E35)-(F38^2)+(2*E35*F38)+(E38*H38))/(K34))+((F42)/(F35*(2-2*G35)));	
N34=((((2*F35+2*E35-G38-H38)*(-M34^2+L34^2))/(J34))+((2*(M34-L34)*(E35-F38)*(E35-F38))/(J34^2)))^(-1);	
O34=(((2*F35)+(2*E35)-G38-H38)*L34/(J34))-(((E35-F38)*(E35-F38))/(J34^2));	
P34=(((2*F35)+(2*E35)-G38-H38)*M34/(J34))-(((E35-F38)*(E35-F38))/(J34^2));	
%
H42= (E35+2*F35)-(F15*(-M34^2+L34^2)*N34);	
I42=(E35)+(F15*(((E35-F38)*(M34-L34))/(J34))*N34);	
J42=(E35+2*F35)-(F15*O34*N34);	
K42=(E35)+(F15*P34*N34);	
L42=(F35)-F15*(((2)/((2*F35)-G38+H38))-(((2*D42)-((4*F42)/(2-2*G35)))/(F35)))^(-1);	
M42=(F35)-F15*((1/(F35-I38))-(D42/F35))^(-1);
% C*
H47=H42;	
I47=I42;	
J47=(0.75*J42+0.25*K42+0.5*M42);	
K47=(0.25*J42+0.75*K42-0.5*M42);	
L47=0.5*(J42-K42);	
M47=M42;
% calculating engineering properties
B24=(H47)-((2*I47*I47)/(J47+K47));	
C24=(2*H47*J47+2*H47*K47-4*I47*I47)*(J47-K47+2*L47)/(3*H47*J47+H47*K47+2*H47*L47-4*I47*I47);	
D24=I47/(J47+K47);	
F24=M47;	
G24=0.25*J47-0.25*K47+0.5*L47;	
E24=C24/(2*G24)-1;	
H24=1-(D24*D24*C24/B24);	
I24=1-(D24*D24*C24/B24);	
J24=1-(E24*E24);	
K24=1-(D24*D24*C24/B24)-(E24*E24)-(D24*D24*C24/B24)-(2*(D24*C24/B24)*(E24)*(D24));
% Engineering properties
E_1=B24;
E_2=C24;
n_12=D24;
n_23=E24;
G_12=F24;
G_23=G24;
if ((H12<=0) | (I12<=0) | (J12<=0) | (K12<=0))
    disp('Error: restriction on OUTPUT elastic constants')
end

% disp('[E_1,E_2,n_12,n_23,G_12,G_23]'); disp([E_1,E_2,n_12,n_23,G_12,G_23]);
fprintf(' E_1 = %g\n E_2 = %g\n n_12 = %g\n n_23 = %g\n G_12 = %g\n G_23 = %g\n',E_1,E_2,n_12,n_23,G_12,G_23)
% end