% FEAcomp. Solution of Example 8.3.
%
% J.A. Mayugo, WVU, 2006

clear all; close all;format short g;

% Volume fraction 
V_f = 0.6;

%% Matrix Epoxy 
E_m=4.5e3; % MPa  

%% Fiber Carbon: 
E_f	= 230e3; %MPa  
m = 8.9;

F1t = 1860 ; % MPa

D_cr = 1- exp(-1/m);

sigma_cr = F1t / (V_f + E_m/E_f *(1-V_f))

delta_alpha = sigma_cr^(-m) / (m*exp(1));
alpha_L = (gamma(1+1/m)/3.53e3)^(m);
delta_Linv = delta_alpha/alpha_L;


disp('Fiber properties');
    disp('E, m, delta_alpha');
    disp([E_f  m  delta_alpha]);
    tline5 = ['!        E      m      delta_alpha'];
    tline7 = ['TBDATA,,' num2str(E_f) ',' num2str(m) ',' num2str(delta_alpha)]; 

n_file = 'FEAcomp_Ex803_to_ANSYS';

tline1 = ['TB,USER,1,1,3, ! Material properties #1, 3 variables'];
tline2 = ['TBTEMP,0'];
tline3 = ['!    Variable descriptions fot the USERMAT sobroutine'];
tline4 = ['!    ------------------------------------------------'];

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

fprintf(fid0,'%10s\n\n','INPUT FILE TO ANSYS, one-dimensional damage using USERMAT');

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',tline7);

fclose(fid0);

% FEA Graph solution
% time, UX, FX
FEsolu =    [   0.20000E-01    0.500000E-02   115.000  ;  
  0.40000E-01    0.100000E-01   230.000 ;   
  0.60000E-01    0.150000E-01   345.000 ;   
  0.80000E-01    0.200000E-01   460.000 ;   
  0.10000        0.250000E-01   575.000 ;   
  0.12000        0.300000E-01   690.000 ;   
  0.14000        0.350000E-01   805.000 ;   
  0.16000        0.400000E-01   919.999 ;   
  0.18000        0.450000E-01   1035.00 ;   
  0.20000        0.500000E-01   1149.99 ;   
  0.22000        0.550000E-01   1264.98 ;   
  0.24000        0.600000E-01   1379.95 ;   
  0.26000        0.650000E-01   1494.89 ;   
  0.28000        0.700000E-01   1609.78 ;   
  0.30000        0.750000E-01   1724.57 ;   
  0.32000        0.800000E-01   1839.18 ;   
  0.34000        0.850000E-01   1953.50 ;   
  0.36000        0.900000E-01   2067.36 ;   
  0.38000        0.950000E-01   2180.50 ;   
  0.40000        0.100000       2292.52 ;   
  0.42000        0.105000       2402.89 ;   
  0.44000        0.110000       2510.83 ;   
  0.46000        0.115000       2615.29 ;   
  0.48000        0.120000       2714.84 ;   
  0.50000        0.125000       2807.60 ;   
  0.52000        0.130000       2891.11 ;   
  0.54000        0.135000       2962.27 ;   
  0.56000        0.140000       3017.23 ;   
  0.58000        0.145000       3051.35 ;   
  0.60000        0.150000       3059.29 ;   
  0.62000        0.155000       3035.09 ;   
  0.64000        0.160000       2972.62 ;   
  0.66000        0.165000       2866.14 ;   
  0.68000        0.170000       2711.19 ;   
  0.70000        0.175000       2505.80 ;   
  0.72000        0.180000       2251.81 ;   
  0.74000        0.185000       1956.13  ;  
  0.76000        0.190000       1631.43  ;  
  0.78000        0.195000       1295.85  ;  
  0.80000        0.200000       971.122  ;  
  0.82000        0.205000       679.165  ;  
  0.84000        0.210000       437.673  ;  
  0.86000        0.215000       256.115  ;  
  0.88000        0.220000       133.814  ;  
  0.90000        0.225000       61.2281  ;  
  0.92000        0.230000       23.9975  ;  
  0.94000        0.235000       7.85579  ;  
  0.96000        0.240000       2.08718  ;  
  0.98000        0.245000      0.435660   ; 
   1.0000        0.250000      0.688618E-01];

[m,n]=size (FEsolu');

time  = FEsolu(:,1)';
epsilon = (FEsolu(:,2)/10*100)';
sigma   = FEsolu(:,3)';

slope_cr = (1-D_cr)*E_f*epsilon/100;

szefont=20;
f1=figure;
plot([0 epsilon],[0 sigma],'-','LineWidth',2.4); 
hold on;
plot([0 epsilon(1:30)],[0 slope_cr(1:30)],'--','LineWidth',1); 
plot(1.5,sigma_cr,'o','LineWidth',1.5,'MarkerSize',10);
%legend('Laminate','0 degree lamina','90 degree lamina','Location','SouthWest');legend('boxoff');
xlabel('\epsilon [%]','Fontsize',szefont);ylabel('\sigma [MPa]','Fontsize',szefont);
linetxt1 = ['    \leftarrow \sigma_{cr}=3060 MPa']
linetxt2 = ['\angle (1-D_{cr})E_f']
linetxt3 = ['    where D_{cr}=10.63%']
text(1.5,sigma_cr,linetxt1,'HorizontalAlignment','left','Fontsize',szefont)
text(0.8,1500,linetxt2,'HorizontalAlignment','left','Fontsize',szefont)
text(0.8,1200,linetxt3,'HorizontalAlignment','left','Fontsize',szefont)
set(gca,'Fontsize',szefont);
hold off;
saveas(f1,['figEx83'],'eps');

