MATLAB Examples

Compression Testing Analysis

compression of a cylindrical gel with height = L [mm] and diameter = d [mm] by Roche C. de Guzman, Ph.D. Hofstra University

Contents

Clear Previous

clear; clc; close('all');

Given

Data = xlsread('DataDispLoad'); % read Excel file: DataDispLoad.xlsx
dL = Data(:,1)'; % relative displacement [mm] data
F = Data(:,2)'; % load = force [N] data
ADF = 10; % adjustment factor for zeroing
L = 1.5875; % height of the cylinder [mm]
d = 4.7625; % diameter of the cylinder [mm]

Computations

% get the zero point of force
NF = numel(F); % number of elements of F
PosF_L = F > 0; % logical positive
IndV = find(PosF_L); % index of positives
NIndV = numel(IndV); % number of elements of IndV
Nneg = NF-NIndV; % number of elements of negatives
Diff = zeros(1,NIndV-ADF);
for c1 = 1:NIndV-ADF
    Diff(c1) = IndV(c1+ADF)-IndV(c1); % difference between points separated by ADF
end
Diff2 = Diff(2:end)-Diff(1:end-1); % difference between 2 adjacent points
for c2 = 1:c1-1
    if Diff2(c2) ~= 0
        Diff2(c2) = c2; % replaces the elements with increasing counter
    end
end
IndZero = max(Diff2)+Nneg+1; % index of zero F

% zeroed dl and F
dL = dL(IndZero:end)-dL(IndZero); % relative displacement [mm]
F = F(IndZero:end)-F(IndZero); % force [N]

% intensive properties
r = d/2; % radius [mm]
A = pi*r^2; % cross-sectional area [mm^2]
A = A/1e6; % cross-sectional area [m^2]
strain = dL/L; % strain [mm/mm]
stress = F/A; % stress [N/m^2 = Pa]

% adjust units
strain = strain*100; % strain [%]
stress = stress/1e6; % stress [MPa]

% locate the ultimate stress
dsigma = diff(stress); % derivative of stress
dsigma = round(dsigma*1000)/1000; % round off dsigma
Ndsig = numel(dsigma); % number of elements of dsigma
MSp = zeros(1,Ndsig); MSpi = MSp; UsigInd = Ndsig; % initial values
for c3 = 1:Ndsig
    [MSp(c3),MSpi(c3)] = max(stress(1:c3)); % maximum stress and its index vectors
    if c3 > 2
        % condition: max is established when the next two are decreasing
        % and the derivative must be negative
        if (MSp(c3) == MSp(c3-1)) && (MSp(c3) == MSp(c3-2)) && (dsigma(c3-1) < 0) && (dsigma(c3-2) < 0)
            UsigInd = c3-2; % index of maximum point
            break;
        end
    end
end

% trim data
cutoff = round(UsigInd*1.15); % maximum index
if cutoff <= Ndsig
    strain = strain(1:cutoff); % strain [%]
    stress = stress(1:cutoff); % stress [MPa]
end

% readjust the zero
Ustr = max(stress); % maximum stress [MPa]
limStr = Ustr*0.01; % 1% of maximum stress [MPa]
Lsig = stress >= limStr; % logical true
LsigInd = find(Lsig); % index of true
zInd = LsigInd(1); % new index
strain = strain(zInd:end)-strain(zInd); % strain [%]
stress = stress(zInd:end)-stress(zInd); % stress [MPa]

% true values
strainT = -log(1-strain/100)*100; % true strain [%]
stressT = stress.*(1-strain/100); % true stress [MPa]
[UCS,UCSi] = max(stressT); % ultimate compressive strength [MPa]
UCstrain = strainT(UCSi); % ultimate compressive strain [%]

% linear regression to determine the compressive modulus
ND = numel(stressT); % number of elements
% initial values
RSQ = zeros(1,ND-2); x = NaN(1,ND); y = x; yf = y; m = RSQ; b = m; dy = b;
for c4 = 1:ND-2
    x = strainT(1:c4+2); % x observed [%]
    y = stressT(1:c4+2); % y observed [MPa]
    m(c4) = (((c4+2)*sum(x.*y))-(sum(x)*sum(y)))/(((c4+2)*sum(x.^2))-(sum(x)^2)) ; % slope of line fit [MPa/% = 100*MPa]
    b(c4) = ((sum(y))-(m(c4)*(sum(x))))/(c4+2); % y-intercept of line fit [MPa]
    yf = m(c4)*x + b(c4); % y fit [MPa]
    SSE = sum((y-yf).^2); % sum of squares error
    SST = sum((y-mean(y)).^2); % sum of squares total
    RSQ(c4) = 1 - (SSE/SST); % coefficient of determination vector
    %xi(c4) = -b/m(c4);
    dy(c4) = abs(y(c4+2) - (m(c4)*x(c4+2)+b(c4)))/UCS; % change in y over UCS
end
[~,ECi] = max(dy >= 0.01); % index of the compressive modulus, 10% cutoff
rsq = RSQ(ECi); % r^2 = coefficient of determination
EC = m(ECi)*100; % compressive modulus [MPa]

% yield values
Ystrain = strainT(ECi+2); % yield strain [%]
Ystress = stressT(ECi+2); % yield stress [%]

Display Results

% animation to fit the line
for c5 = 1:ND-2
    plot(strain,stress,'-r'); % engineering
    hold('on');
    plot(strainT,stressT,'-b','linewidth',2); % true
    legend('engineering','true','location','northwest');
    title('Gel Compression Until Fracture');
    xlabel('Compressive Strain [%]');
    ylabel('Compressive Stress [MPa]');
    plot([strainT(1) strainT(end)],[m(c5)*strainT(1)+b(c5) m(c5)*strainT(end)+b(c5)],'color',[0 1 0.25]);
    axis([0 strainT(ND-2) 0 max(stress)*1.05]);
    drawnow;
    hold('off');
end
for c6 = ND-2:-1:1
    plot(strain,stress,'-r'); % engineering
    hold('on');
    plot(strainT,stressT,'-b','linewidth',2); % true
    legend('engineering','true','location','northwest');
    title('Gel Compression Until Fracture');
    xlabel('Compressive Strain [%]');
    ylabel('Compressive Stress [MPa]');
    plot([strainT(1) strainT(end)],[m(c6)*strainT(1)+b(c6) m(c6)*strainT(end)+b(c6)],'color',[0 1 0.25]);
    axis([0 strainT(ND-2) 0 max(stress)*1.05]);
    drawnow;
    hold('off');
    if m(c6)*100 == EC
        plot(strain,stress,'-r'); % engineering
        hold('on');
        plot(strainT,stressT,'-b','linewidth',2); % true
        legend('engineering','true','location','northwest');
        title('Gel Compression Until Fracture');
        xlabel('Compressive Strain [%]');
        ylabel('Compressive Stress [MPa]');
        % yield
        plot([strainT(1) strainT(end)],[m(c6)*strainT(1)+b(c6) m(c6)*strainT(end)+b(c6)],'-g','color',[0 1 0.25]);
        plot(Ystrain,Ystress,'ob','markerfacecolor',[0 1 0],'markersize',5);
        plot([Ystrain Ystrain],[0 Ystress],'--g');
        plot([0 Ystrain],[Ystress Ystress],'--g');
        text(Ystrain*0.8,Ystress+0.08*UCS,'yield','color',[0 1 0]);
        % ultimate
        plot(UCstrain,UCS,'pb','markerfacecolor',[0 0 0],'markersize',8);
        plot([UCstrain UCstrain],[0 UCS],'--k');
        plot([0 UCstrain],[UCS UCS],'--k');
        text(UCstrain*1.01,UCS+0.08*UCS,'ultimate');
        axis([0 strainT(ND-2) 0 max(stress)*1.05]);
        hold('off');
        break;
    end
end

% command window display
disp('========================================================================================');
disp(['   The gel was found to have a compressive modulus of ' num2str(EC) ' MPa or ' num2str(EC*1000) ' kPa.']);
disp(['   Its yield strain is ' num2str(Ystrain) '%, while']);
disp(['   its yield strength is ' num2str(Ystress) ' MPa or ' num2str(Ystress*1000) ' kPa.']);
disp(['   Also, its ultimate compressive strain is ' num2str(UCstrain) '%, while']);
disp(['   its ultimate compressive strength is ' num2str(UCS) ' MPa or ' num2str(UCS*1000) ' kPa.']);
disp('========================================================================================');

% write to Excel
Results = {'Compressive Modulus of Elasticity [MPa]' 'Yield Strain [%]' 'Yield Stress [MPa]' 'Ultimate Compressive Strain [%]' 'Ultimate Compressive Stress [MPa]'; EC Ystrain Ystress UCstrain UCS}; % summary of results
xlswrite('ResStressStrain.xlsx',Results); % write to Excel file: ResStressStrain.xlsx

% Thank You!!!
========================================================================================
   The gel was found to have a compressive modulus of 0.1836 MPa or 183.5993 kPa.
   Its yield strain is 8.5438%, while
   its yield strength is 0.015878 MPa or 15.8781 kPa.
   Also, its ultimate compressive strain is 37.8399%, while
   its ultimate compressive strength is 0.16715 MPa or 167.1477 kPa.
========================================================================================