% FEAcomp. Solution of Example 7.4

% To run this M-file need the functions:
%       fitfunFEAcomp(xdata, ydata,ifun)
%       invlapFEAcomp(time,func)
%       PMMViscoMatrix( time, V_f, imod_m, imod_f, Elas_prop_m, Elas_prop_f, Visc_prop_m)
% available in .\Examples\Ch7ex\FEAcomp_Ex704\

clear all; close all;format short g;

%% domain time
time = [0.001;.01;.1;1;2;3;5;10;20;30;50;100;200;300];

% Volume fraction 
V_f = 0.6;

%% Matrix Vinylester 
E_m=4.08; % GPa 0.592;      % msi 
nu_m = 0.311;

% Maxwell model
imod_m=1;
eta_m = 159.8; %GPa/min23.18;  %msi/min

%%% Fiber: 
for imod_f=2:2;  % 1,1 only isotropic, 2:2 only transversaly, 1:2 both
    if imod_f==1 % isotropic
        E_f	= 168.4; %GPa =24.42 msi 
        nu_f  = 0.443;	
    elseif imod_f==2 % transversely isotropic																					
        E_a  = 168.4;  E_t = 24.8; % GPa = 3.6 msi 
        nu_a = 0.443;  nu_t = 0.005; 
        G_a  = 44.1; % GPa =6.4 msi;    
        G_t = E_t/2/(1+nu_t); % GPa, calculated
    end

    if imod_f==1       % Isotropic fibers, requires imod_f =1
        Elas_prop_f=[E_f nu_f];     
    elseif imod_f==2    % Transversaly Isotropic fibers, requires imod_f =2
        Elas_prop_f=[E_a E_t nu_a nu_t G_a G_t];     
    end
    if imod_m==1         % Maxell model, requires imod_m = 1
        Visc_prop_m=[eta_m];
    else
        Visc_prop_m=[0];
    end
    Elas_prop_m=[E_m nu_m];
    [Elas_prop, Visc_prop,model] = PMMViscoMatrix( time, V_f, imod_m, imod_f,...
                                             Elas_prop_m, Elas_prop_f, Visc_prop_m);
    %% Elastic properties
    Eo_1=Elas_prop(1); Eo_2=Elas_prop(2);
    Eo_12=Elas_prop(3); Eo_23=Elas_prop(4);
    nuo_12=Elas_prop(5); nuo_23=Elas_prop(6);

    disp('Composite: Elastic properties');
    disp('Eo_1, Eo_2, Go_12, Go_23, nuo_12, nuo_23');
    disp([Eo_1 Eo_2 Eo_12 Eo_23 nuo_12 nuo_23]);
    tline5 = ['!        E1     E2     nu12    nu23    G12     G23'];
    tline7 = ['TBDATA,,' num2str(Eo_1) ',' num2str(Eo_2) ',' num2str(nuo_12) ',' num2str(nuo_23) ',' num2str(Eo_12) ',' num2str(Eo_23)];
    
    %% Viscoelastic properties                                     
    if imod_m~=0
        if imod_m==1         % Maxell model, requires imod_m = 1
            tau_1=Visc_prop(1); tau_2=Visc_prop(2); tau_12=Visc_prop(3); tau_23=Visc_prop(4);         
        
        %% Create exponential fit and plot E_2t in a separate figure
        figure1=figure;
        %[sse, FittedCurve] = model([Eo_2 1/tau_2]);
        %plot((time), FittedCurve,'x','MarkerSize',16,'LineWidth',1.5);
        hold on;
        dtime=[0.001:100];
        plot((dtime),Eo_2.*exp(-dtime./tau_2),'k','LineWidth',2.4)
        plot((dtime),Eo_2.*(dtime./dtime),'k--','LineWidth',1.2)
        xlabel('t [minutes]','Fontsize',16);ylabel('E_2(t) [GPa]','Fontsize',16);
        legend('Viscoelastic response','Elastic','Location','SouthWest');legend('boxoff');
        axis([0 100 0 15])
        set(gca,'Fontsize',16);
        
        saveas(figure1,['Fig7_4'],'eps');

        disp('Composite: relaxations as follows (Poison''s remain constant)');
        disp('tau_1, tau_2, tau_12, tau_23');
        disp([tau_1 tau_2 tau_12 tau_23]);
        
        tline6 = ['!        Tau1   Tau2   Tau12   Tau13'];
        tline8 = ['TBDATA,,' num2str(tau_1) ',' num2str(tau_2) ',' num2str(tau_12) ',' num2str(tau_23)];
        end
    else
        tline6 = [' '];  tline8 = [' '];
    end
end
hold off;

n_file = 'FEAcomp_Ex704_to_ANSYS';

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

% disp(tline1);
% disp(tline2);
% disp(tline3);
% disp(tline4);
% disp(tline5);
% disp(tline6);
% disp(tline7);
% disp(tline8);

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

fprintf(fid0,'%10s\n\n','INPUT FILE TO ANSYS, viscoelasticity 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',tline6);
fprintf(fid0,'%10s\n',tline7);
fprintf(fid0,'%10s\n',tline8);

fclose(fid0);

