how can I calculate rate of compression?

9 views (last 30 days)
Manav Divekar
Manav Divekar on 30 Jan 2022
Answered: Vinayak on 11 Jan 2024
I want to test rate of compression (mm/s), average Elastic moduli and Plot force-displacement. from the given given text file. Since there are other text files whoes block size are different what changes can i make in my code to make it work for all without manually changing values. For this particula set of value 10% strain rate was applied.
Data = inportdata('DHD_10%_sec_T1.txt'); % read text file: DHD_10%_sec_T1.txt
dL = Data(:,1)'; % relative displacement [mm] data
F = Data(:,2)'; % load = force [N] data
ADF = 10; % adjustment factor for zeroing
L = 39.96; % lenght if block [mm]
w = 25.8666667; % width if block[mm]
th = 25.9; % thickness if block[mm]
% 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
A = w*L*th; % area [mm^2]
A = A/1e6; % 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 [%]

Answers (1)

Vinayak
Vinayak on 11 Jan 2024
Hi Manav
From what I noticed in the file you shared, it contains three columns of data. I presume that the other files are in similar format. You could write a custom function that determines the value of strain rate and block sizes based on file name. You can also write a function to go through all files in an existing directory and process them accordingly. Here is a sample code for the given file to calculate the strain rate:
% Get a list of all files in the directory
directoryPath = 'path/to/your/directory'; % Replace with your directory path
files = dir(directoryPath);
fileNames = {files.name};
for i = 1:length(fileNames)
fileName = fileNames{i};
% Ignore directories and hidden files
if ~files(i).isdir && ~startsWith(fileName, '.')
% Extract the strain rate from the file name
strainRate = getStrainRateFromFileName(fileName);
% Display the strain rate
fprintf('Processing file: %s with strain rate: %d%%\n', fileName, strainRate);
% Process the stress/strain data from the file
processData(fullfile(directoryPath, fileName), strainRate);
end
end
function strainRate = getStrainRateFromFileName(fileName)
% Regex pattern to match digits followed by a percentage sign
pattern = '_(\d+)%';
tokens = regexp(fileName, pattern, 'tokens');
if ~isempty(tokens)
strainRate = str2double(tokens{1}{1});
else
error('Strain rate could not be extracted from the file name.');
end
end
function processData(filePath, strainRate)
Data = importdata(filePath);
ADF = strainRate;
% ... (rest of the old code)
end
You may read more about regexp and dir” in the following documentation:

Categories

Find more on Stress and Strain in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!