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

lm_0=E_m*nu_m/((1+nu_m)*(1-2*nu_m));	
mu_0=E_m/(2*(1+nu_m));	
nu_0=nu_m;		

%%%Fiber: isotropic																					
E_f	= 24.42; nu_f  = 0.443;	

lm_1=nu_f*E_f/((1+nu_f)*(1-2*nu_f));	
mu_1=E_f/(2*(1+nu_f));	
nu_1=nu_f;																					

%% a	b	c	g

a=mu_1-mu_0-2*mu_1*nu_m+2*mu_0*nu_f;	
b=-mu_0*nu_m+mu_1*nu_f+2*mu_0*nu_m*nu_f-2*mu_1*nu_m*nu_f;	
c=(mu_0-mu_1)*(-mu_0+mu_1-mu_0*nu_m-2*mu_1*nu_m+2*mu_0*nu_f+mu_1*nu_f+2*mu_0*nu_m*nu_f-2*mu_1*nu_m*nu_f);	
g=2-2*nu_m;
	
V_f = 0.6;	

%%S_3	S_6	S_7
S_3=0.49247-0.47603*V_f-0.02748*V_f^2;	
S_6=0.36844-0.14944*V_f-0.27152*V_f^2;	
S_7=0.12346-0.32035*V_f+0.23517*V_f^2;	

%%D_	C_11*	C_12*	C_23*	C_22*	C_44*	C_66*
D_=(a*S_3^2)/(2*mu_0^2*c)-(a*S_6*S_3)/(mu_0^2*g*c)+a*(S_6^2-S_7^2)/(2*mu_0^2*g^2*c)+S_3*(b^2-a^2)/(2*mu_0*c^2)+(S_6*(a^2-b^2)+S_7*(a*b+b^2))/(2*mu_0*g*c^2)+(a^3-2*b^3-3*a*b^2)/(8*c^3);	
C_11=lm_0+2*mu_0-V_f*((S_3^2)/(mu_0^2)-2*S_6*S_3/(mu_0^2*g)-a*S_3/(mu_0*c)+(S_6^2-S_7^2)/(mu_0^2*g^2)+(a*S_6+S_7*b)/(mu_0*g*c)+(a^2-b^2)/(4*c^2))/D_;	
C_12=lm_0+((S_3/(2*c*mu_0))-((S_6-S_7)/(2*mu_0*c*g))-(a+b)/(4*c^2))*b*V_f/D_;	
C_23=lm_0+V_f*((a*S_7)/(2*mu_0*g*c)-(b*a+b^2)/(4*c^2))/D_;	
C_22=lm_0+2*mu_0-V_f*((-a*S_3)/(2*mu_0*c)+(a*S_6)/(2*mu_0*g*c)+(a^2-b^2)/(4*c^2))/D_;	
C_44=mu_0-V_f/((-2*S_3)/mu_0+1/(mu_0-mu_1)+4*S_7/(mu_0*(2-2*nu_0)));	
C_66=mu_0-V_f/((-S_3/mu_0)+1/(mu_0-mu_1));	

%%% Engineering Constants
%%E_1	E_2	n_12	n_23	G_12	G_23		
E_1=C_11-(2*C_12^2)/(C_22+C_23);	
E_2=((2*C_11*C_22+2*C_11*C_23-4*C_12^2)*(C_22-C_23+2*C_44))/(3*C_11*C_22+C_11*C_23+2*C_11*C_44-4*C_12^2);	
n_12=C_12/(C_22+C_23);	
n_23=(C_11*C_22+3*C_11*C_23-2*C_11*C_44-4*C_12^2)/(3*C_11*C_22+C_11*C_23+2*C_11*C_44-4*C_12^2);	
G_12=C_66;	
G_23=E_2/(2*(1+n_23));

% 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)

%%% Restrictions on elastic constants
%%1-n12n21>0	1-n13n31>0	1-n23n32>0	1-n12n21-n23n32-n13n31-2n21n32n13>0
%=1-(n_12.*n_12.*E_2./E_1)	
%=1-(n_12.*n_12.*E_2./E_1)	
%=1-(n_23.*n_23)	
%=1-(n_12.*n_12.*E_2./E_1)-(n_23.*n_23)-(n_12.*n_12.*E_2./E_1)-(2.*(n_12.*E_2./E_1).*(n_23).*(n_12))					

% end
