diff --git a/Fingertip.m b/Fingertip.m
index b626037dcc44b0dd2a23f2c242d469e932e8e19a..1eea1c1dc7b6e5dca5d7a9aec4157e078241f189 100644
--- a/Fingertip.m
+++ b/Fingertip.m
@@ -16,7 +16,7 @@ classdef Fingertip
b = 0.1;
% % Sampling frequency
- Fs = 4e5;
+ Fs = 6e5;
% % External forces applied on the bone
P = 2; %% initial pressure
@@ -34,7 +34,7 @@ classdef Fingertip
end
methods
function Ftip = Init(Ftip,varargin)
- global kext kt kn Kn B x0
+ global kext kt kn Kn B x0 l0
N = Ftip.Nl;
@@ -50,7 +50,7 @@ classdef Fingertip
Ftip.P = p.Results.pressure;
Ftip.Q = p.Results.traction;
- Ftip.kc = p.Results.stiffness*5e-3;
+ Ftip.kc = p.Results.stiffness*Ftip.d/2;
Ftip.fclb = p.Results.friction;
% Ftip.nu = min(p.Results.friction,0.48);
Ftip.Rcurv = p.Results.curvature;
@@ -78,7 +78,7 @@ classdef Fingertip
Kn = sparse(idx_i,idx_j,values,2*(N+1),2*(N+1));
end
- function dU = f_U(Ftip,Uc,Fc)
+ function [dU,SE] = f_U(Ftip,Uc,Fc)
global kext kt kn Kn B x0
N = Ftip.Nl;
pos_x = Uc(N+2:2*(N+1))+x0(N+2:2*(N+1));
@@ -107,8 +107,8 @@ classdef Fingertip
d21sub = [0,part_ext.*Asin(2:N),0]';
d21sup = [0,0,part_ext.*Asin(2:N)]';
spKd = (spdiags([d0sub d0 d0sup],-1:1,length(d0),length(d0)))+...
- 0*(spdiags([zeros(N+1,3);[d12sub d12 d12sup]],N:N+2,length(d0),length(d0)))+...
- 0*(spdiags([[d21sub d21 d21sup];zeros(N+1,3)],-N-2:-N,length(d0),length(d0)));
+ Ftip.nu/100*(spdiags([zeros(N+1,3);[d12sub d12 d12sup]],N:N+2,length(d0),length(d0)))+...
+ Ftip.nu/100*(spdiags([[d21sub d21 d21sup];zeros(N+1,3)],-N-2:-N,length(d0),length(d0)));
% Ksh: Springs dependencies of the subcutaneous tissues
angleB = atan(abs(pos_z(1)-pos_z(2:N+1))./abs(pos_x(1)-pos_x(2:N+1))); %% angle between two elements linked to the bone
@@ -136,10 +136,12 @@ classdef Fingertip
spK = kext.*spKd+kt.*spKsh+kn.*Kn;
dU = -inv(B)*spK*Uc-inv(B)*Fc;
+
+ SE = 1/2*spK*Uc.^2;
end
- function [U,F] = Run(Ftip)
- global x0
+ function [U,F,SE,eps] = Run(Ftip,deep)
+ global x0 l0
%% Time vector
dt = 1/Ftip.Fs; %% sampling period
time = 0:dt:0.5; %% time vector
@@ -160,6 +162,9 @@ classdef Fingertip
Upt = zeros(N+1,1); %% speed tangential to the surface
pos = x0; %% position of all elements
+ ndraw = 500;
+ X = x0(N+3):1e-4:x0(end); Zv = (-4e-3:2e-4:-1e-4)';
+ [Xq,Zq]=meshgrid(X,Zv);
tic
while(t<=2 || abs(sum(F(1:N+1,t-1),1))>10e-2*Ftip.P) %% Stop condition: global net force <5%
if(t>2)
@@ -186,13 +191,53 @@ classdef Fingertip
k1 = f_U(Ftip,U(:,t-1),F(:,t));
k2 = f_U(Ftip,U(:,t-1)+0.5*dt*k1,F(:,t));
k3 = f_U(Ftip,U(:,t-1)+0.5*dt*k2,F(:,t));
- k4 = f_U(Ftip,U(:,t-1)+dt*k3,F(:,t));
+ [k4,SE_t] = f_U(Ftip,U(:,t-1)+dt*k3,F(:,t));
U(:,t) = U(:,t-1)+ dt/6*(k1+2*k2+2*k3+k4);
-
+ SE(:,t) = SE_t; %mJ/mm^3
pos = U(:,t) + x0;
+ if(deep == true)
+ if(~mod(t,ndraw))
+ %% Strains at the depth of the mechanoreceptors
+ p(:,t-1) = F(2:N+1,t-1)./(l0*Ftip.d);
+ q(:,t-1) = F(N+3:2*(N+1),t-1)./(l0*Ftip.d);
+ pf = @(s) p(find(abs(x0(N+3:2*(N+1))-s)==min(abs(x0(N+3:2*(N+1))-s))),t-1);
+ qf = @(s) q(find(abs(x0(N+3:2*(N+1))-s)==min(abs(x0(N+3:2*(N+1))-s))),t-1);
+
+ sigma_x = -2/pi*Zq.*integral(@(s) pf(s).*(Xq-s).^2./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true)-...
+ -2/pi.*integral(@(s) qf(s).*(Xq-s).^3./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true);
+ sigma_z = -2/pi*Zq.^3.*integral(@(s) pf(s)./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true)-...
+ -2/pi*Zq.^2.*integral(@(s) qf(s).*(Xq-s)./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true);
+% tau = @(x) -2/pi*Zq.^2.*integral(@(s) pf(s).*(Xq-s).^2./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true)-...
+% -2/pi*Zq.*integral(@(s) qf(s).*(Xq-s).^2./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true);
+ eps(:,:,t) = 1/Ftip.E_tis*[sigma_x-Ftip.nu*sigma_z;-Ftip.nu*sigma_x+sigma_z];
+ end
+ else
+ eps(:,:,t) = zeros(size(Xq));
+ end
t = t+1;
end
+ toc
+ if(deep == true)
+ %% Strains at the depth of the mechanoreceptors
+ tic
+ p(:,t-1) = F(2:N+1,t-1)./(l0*Ftip.d);
+ q(:,t-1) = F(N+3:2*(N+1),t-1)./(l0*Ftip.d);
+ pf = @(s) p(find(abs(x0(N+3:2*(N+1))-s)==min(abs(x0(N+3:2*(N+1))-s))),t-1);
+ qf = @(s) q(find(abs(x0(N+3:2*(N+1))-s)==min(abs(x0(N+3:2*(N+1))-s))),t-1);
+
+ sigma_x = -2/pi*Zq.*integral(@(s) pf(s).*(Xq-s).^2./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true)-...
+ -2/pi.*integral(@(s) qf(s).*(Xq-s).^3./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true);
+ sigma_z = -2/pi*Zq.^3.*integral(@(s) pf(s)./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true)-...
+ -2/pi*Zq.^2.*integral(@(s) qf(s).*(Xq-s)./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true);
+% tau = @(x) -2/pi*Zq.^2.*integral(@(s) pf(s).*(Xq-s).^2./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true)-...
+% -2/pi*Zq.*integral(@(s) qf(s).*(Xq-s).^2./((Xq-s).^2+Zq.^2).^2,x0(N+3),x0(end),'ArrayValued',true);
+ eps(:,:,t-1) = 1/Ftip.E_tis*[sigma_x-Ftip.nu*sigma_z;-Ftip.nu*sigma_x+sigma_z];
+ toc
+ else
+ eps(:,:,t) = zeros(size(Xq));
+ end
+
ti = t;
if(Ftip.Q~=0)
F(:,ti:Ndt) = repmat([-Ftip.P;zeros(N,1);Ftip.Q;zeros(N,1)],1,Ndt-ti+1);
@@ -219,32 +264,34 @@ classdef Fingertip
k1 = f_U(Ftip,U(:,t-1),F(:,t));
k2 = f_U(Ftip,U(:,t-1)+0.5*dt*k1,F(:,t));
k3 = f_U(Ftip,U(:,t-1)+0.5*dt*k2,F(:,t));
- k4 = f_U(Ftip,U(:,t-1)+dt*k3,F(:,t));
+ [k4,SE_t] = f_U(Ftip,U(:,t-1)+dt*k3,F(:,t));
U(:,t) = U(:,t-1)+ dt/6*(k1+2*k2+2*k3+k4);
-
+ SE(:,t) = SE_t;
pos = U(:,t) + x0;
t = t+1;
end
end
- toc
U = U(:,1:t-1);
F = F(:,1:t-1);
end
- function [] = Animate(Ftip,U,F,save)
+ function [] = Animate(Ftip,U,F,eps,deep,save)
global x0
N = Ftip.Nl;
%% Surface
surface = @(x) (sign(Ftip.Rcurv)*(sqrt(abs(Ftip.Rcurv)^2-x.^2)-abs(Ftip.Rcurv))).*(Ftip.Rcurv~=0) +...
0.*(Ftip.Rcurv~=0);
- vecSurf = -4e-3:1e-5:4e-3;
+ vecSurf = -5e-3:1e-5:5e-3;
surfDef = surface(vecSurf);
mu = abs(F(N+2:2*(N+1),:)./F(1:N+1,:));
- h = figure('units','normalized','outerposition',[0 0 2/3 1/2],'Color','w');
+ X = x0(N+3):1e-4:x0(end);
+ rho=-3.9e-3:2e-4:0; rho = 0:2e-4:3.9e-3;
+
+ h = figure('units','normalized','outerposition',[0 0 2/3 1],'Color','w');
if(save==true)
filename = ['../GIF/simu_f' num2str(Ftip.fclb) '_c' num2str(Ftip.Rcurv) '_s' num2str(Ftip.kc/5e-3*1e-3) '.gif'];
frame = getframe(h);
@@ -256,14 +303,17 @@ classdef Fingertip
axis off
ndraw = 500;
vec_plot = floor((N+1)/2)+1-69:3:floor((N+1)/2)+1+70;
+ x_m = [-4.92,-3.69,-2.46,-1.23,0,1.23,2.46,3.69,4.92].*1e-3;
sz = 30;
for k=1:size(U,2)
- if (~mod(k-1,ndraw))
+ if (~mod(k,ndraw))
+ if(~deep)
clf;
+ subplot 121
axis off
hold on
% % Update surface deformation
- indC = find(U(1:N+1,k)+x0(1:N+1)<=0);
+ indC = find(U(1:N+1,k)+x0(1:N+1) - surface(U(N+2:2*(N+1),k)+x0(N+2:2*(N+1)))<=0);
if(length(indC)>2)
[~,indS] = min(abs(U(indC+N+1,k)+x0(indC+N+1)-vecSurf)');
surfDef(min(indS):max(indS)) = interp1(indS,U(indC,k)+x0(indC),min(indS):max(indS));
@@ -275,14 +325,97 @@ classdef Fingertip
scatter(U(N+1+vec_plot,k)+x0(N+1+vec_plot),U(vec_plot,k)+x0(vec_plot),sz,cd,'filled');
colorbar;
caxis([0 1])
+ shading flat
+ axis equal
+ axis off
+
+ quiver(U(N+1+vec_plot(2:end),k)+x0(N+1+vec_plot(2:end)),U(vec_plot(2:end),k)+x0(vec_plot(2:end)),...
+ -F(N+1+vec_plot(2:end),k),-F(vec_plot(2:end),k),0.75,'Color','w','Linewidth',1)
+
+ xlim([-5e-3 5e-3])
+ xlabel('x (m)')
+ ylim([-1e-3 3.5e-3])
+ ylabel('y (m)')
+ caxis([-1e-3 3e-3])
+
+ subplot 122
+ hold on
+ pos_x = 1/2*(U(N+1+vec_plot(2:end-1),k)+x0(N+1+vec_plot(2:end-1))+U(N+1+vec_plot(3:end),k)+x0(N+1+vec_plot(3:end)));
+ plot(pos_x.*1e3,diff(U(N+1+vec_plot(2:end),k)))
+ xlabel('position (mm)')
+ ylabel('tangential strain')
+
+
+ else
+ clf;
+ subplot 221
+ axis off
+ hold on
+ % % Update surface deformation
+ indC = find(U(1:N+1,k)+x0(1:N+1) - surface(U(N+2:2*(N+1),k)+x0(N+2:2*(N+1)))<=0);
+ if(length(indC)>2)
+ [~,indS] = min(abs(U(indC+N+1,k)+x0(indC+N+1)-vecSurf)');
+ surfDef(min(indS):max(indS)) = interp1(indS,U(indC,k)+x0(indC),min(indS):max(indS));
+ end
+ area(vecSurf,surfDef,-2e-3,'Facecolor','k')
+
+ z(:,k) = interp1(U(N+3:end,k)+x0(N+3:end),U(2:N+1,k)+x0(2:N+1),X);
+ pcolor(repmat(X,length(rho),1),flipud((z(:,k)+rho)'),eps(1:length(rho),:,k).*1e-3)
+ shading flat
+ axis equal
+ axis off
+
+ z_m = interp1(U(N+3:end,k)+x0(N+3:end),U(2:N+1,k)+x0(2:N+1)+2e-3,x_m); %%PC densitiy=0.81/cm2
+ scatter(x_m,z_m,20,'k','filled')
+
+ quiver(U(N+1+vec_plot(2:end),k)+x0(N+1+vec_plot(2:end)),U(vec_plot(2:end),k)+x0(vec_plot(2:end)),...
+ -F(N+1+vec_plot(2:end),k),-F(vec_plot(2:end),k),0.75,'Color','w','Linewidth',1)
+
+ xlim([-5e-3 5e-3])
+ xlabel('x (m)')
+ ylim([-1e-3 3.5e-3])
+ ylabel('y (m)')
+ caxis([-1e-3 3e-3])
+
+ subplot 222
+ hold on
+ eps_m = interp1(X,eps(find(rho>2e-3,1,'first'),:,k),x_m); %%PC densitiy=0.81/cm2
+ plot(X.*1e3,eps(find(rho>2e-3,1,'first'),:,k))
+ scatter(x_m.*1e3,eps_m,20,'k','filled')
+ ylim([-1 1])
+ xlabel('position (mm)')
+ ylabel('tangential strain')
+
+ subplot 223
+ axis off
+ hold on
+ area(vecSurf,surfDef,-2e-3,'Facecolor','k')
+ z(:,k) = interp1(U(N+3:end,k)+x0(N+3:end),U(2:N+1,k)+x0(2:N+1),X);
+ pcolor(repmat(X,length(rho),1),flipud((z(:,k)+rho)'),eps(1+length(rho):end,:,k).*1e-3)
+ shading flat
+ axis equal
+ axis off
+ z_m = interp1(U(N+3:end,k)+x0(N+3:end),U(2:N+1,k)+x0(2:N+1)+2e-3,x_m); %%PC densitiy=0.81/cm2
+ scatter(x_m,z_m,20,'k','filled')
quiver(U(N+1+vec_plot(2:end),k)+x0(N+1+vec_plot(2:end)),U(vec_plot(2:end),k)+x0(vec_plot(2:end)),...
- F(N+1+vec_plot(2:end),k),F(vec_plot(2:end),k),0.5,'Color','b')
+ -F(N+1+vec_plot(2:end),k),-F(vec_plot(2:end),k),0.75,'Color','w','Linewidth',1)
- xlim([-4e-3 4e-3])
+ xlim([-5e-3 5e-3])
xlabel('x (m)')
- ylim([-1e-3 2e-3])
+ ylim([-1e-3 3.5e-3])
ylabel('y (m)')
+% caxis([-1e-3 3e-3])
+
+ subplot 224
+ hold on
+ eps_m = interp1(X,eps(length(rho)+find(rho>2e-3,1,'first'),:,k),x_m); %%PC densitiy=0.81/cm2
+ plot(X.*1e3,eps(length(rho)+find(rho>2e-3,1,'first'),:,k))
+ scatter(x_m.*1e3,eps_m,20,'k','filled')
+ ylim([-0.5 3])
+ xlabel('position (mm)')
+ ylabel('normal strain')
+ end
drawnow
if(save==true)
@@ -305,9 +438,11 @@ classdef Fingertip
function ac = contactArea(Ftip,U)
global x0
+ surface = @(x) (sign(Ftip.Rcurv)*(sqrt(abs(Ftip.Rcurv)^2-x.^2)-abs(Ftip.Rcurv))).*(Ftip.Rcurv~=0) +...
+ 0.*(Ftip.Rcurv~=0);
N = Ftip.Nl;
for t=1:size(U,2)
- indC = find(U(1:N+1,t)+x0(1:N+1)<=0);
+ indC = find(U(1:N+1,t)+x0(1:N+1) - surface(U(N+2:2*(N+1),t)+x0(N+2:2*(N+1)))<=0);
if(length(indC)>2)
ac(t) = U(max(indC)+N+1,t)+x0(max(indC)+N+1) - (U(min(indC)+N+1,t)+x0(min(indC)+N+1));
end
@@ -321,19 +456,23 @@ classdef Fingertip
function div = divergence(Ftip,U)
global x0
+ surface = @(x) (sign(Ftip.Rcurv)*(sqrt(abs(Ftip.Rcurv)^2-x.^2)-abs(Ftip.Rcurv))).*(Ftip.Rcurv~=0) +...
+ 0.*(Ftip.Rcurv~=0);
N = Ftip.Nl;
- indC = find(U(1:N+1,size(U,2))+x0(1:N+1)<=0);
+ indC = find(U(1:N+1,size(U,2))+x0(1:N+1) - surface(U(N+2:2*(N+1),size(U,2))+x0(N+2:2*(N+1)))<=0);
for t=1:size(U,2)
div2(:,t) = 2*gradient(U(N+3:2*(N+1),t),U(N+3:2*(N+1),t)+x0(N+3:2*(N+1)));
end
- div = nanmedian(div2(indC,:),1);
+ div = nanmedian(div2(:,:),1);
end
function [sedXY,sedXZ] = strainEnergy(Ftip,U)
global x0
+ surface = @(x) (sign(Ftip.Rcurv)*(sqrt(abs(Ftip.Rcurv)^2-x.^2)-abs(Ftip.Rcurv))).*(Ftip.Rcurv~=0) +...
+ 0.*(Ftip.Rcurv~=0);
N = Ftip.Nl;
E = Ftip.E_tis;
- indC = find(U(1:N+1,size(U,2))+x0(1:N+1)<=0);
+ indC = find(U(1:N+1,size(U,2))+x0(1:N+1) - surface(U(N+2:2*(N+1),size(U,2))+x0(N+2:2*(N+1)))<=0);
for t=1:size(U,2)
norm_strain(:,t) = diff(U(2:N+1,t));
tan_strain(:,t) = diff(U(N+3:2*(N+1),t));
@@ -342,8 +481,8 @@ classdef Fingertip
udXZ(:,t) = E*(1-Ftip.nu)/(2*(1+Ftip.nu)*(1-2*Ftip.nu))*(norm_strain(:,t).^2+2*tan_strain(:,t).^2)+...
E*Ftip.nu/((1+Ftip.nu)*(1-2*Ftip.nu))*tan_strain(:,t).^2.*norm_strain(:,t); %linear
- sedXY(:,t) = 0.002*nanmedian(udXY(indC,t)).*1e9; %mJ/mm3
- sedXZ(:,t) = nanmedian(udXZ(indC,t)).*1e9; %mJ/mm3
+ sedXY(:,t) = 0.002*nanmedian(udXY(:,t)).*1e3; %mJ/mm3
+ sedXZ(:,t) = nanmedian(udXZ(:,t)).*1e6; %mJ/mm3
end
end
diff --git a/README.md b/README.md
index 41e76c26dceaaeda4d324cfada83bb0040365230..1950765c6016546250066c50abc62695f2a87ad4 100644
--- a/README.md
+++ b/README.md
@@ -1,7 +1,17 @@
# Finger model
-This code provides a 2D parsimonious finite-difference model of human skin capturing the deformation of the skin, viscoelastic effects through the damper, and local friction behavior at the interface. The model is constructed with a bottom-up approach, using as few parameters to explain a wide array of observed phenomena. It is composed of a chain of massless elements maintained together by elastic springs. This chain can be assimilated to the external layer of the skin (the epidermis). Its shape is maintained using other elastic springs that connect the massless elements to a virtual bone, analogous to the mechanical behavior of the subcutaneous tissues. The two elements on the outside of the membrane are also attached to the bone and model the effect of the rigid nail. Overall, the model resemble the discrete version of a curved elastic membrane on a spring foundation. Viscosity of the skin is modeled by dampers, connecting each particles to the mass of the system.
+This code provides a 2D parsimonious finite-difference model of human skin capturing the deformation of the skin,
+viscoelastic effects through the damper, and local friction behavior at the interface.
+The model is constructed with a bottom-up approach, using as few parameters to explain a wide array of observed phenomena.
+It is composed of a chain of massless elements maintained together by elastic springs.
+This chain can be assimilated to the external layer of the skin (the epidermis).
+Its shape is maintained using other elastic springs that connect the massless elements to a virtual bone,
+analogous to the mechanical behavior of the subcutaneous tissues.
+The two elements on the outside of the membrane are also attached to the bone and model the effect of the rigid nail.
+Overall, the model resemble the discrete version of a curved elastic membrane on a spring foundation.
+Viscosity of the skin is modeled by dampers, connecting each particles to the mass of the system.
-External forces are applied on the bone element and the frictional forces are computed with the Dahl's friction model. The displacements were solved using Runge-Kutta at order 4. For the equations, please refer to the submission.
+External forces are applied on the bone element and the frictional forces are computed with the Dahl's friction model.
+The displacements were solved using Runge-Kutta at order 4. For the equations, please refer to the submission.
diff --git a/example.m b/example.m
index 0e7d317691ae2b71408743b3bd5385728b491815..65514cf8b23e53824f09914357f4ff54cdfa0c7b 100644
--- a/example.m
+++ b/example.m
@@ -2,25 +2,28 @@ clear all
close all
Ftip = Fingertip;
-Ftip = Init(Ftip,'pressure',3);
-[U,F] = Run(Ftip);
-Animate(Ftip,U,F,true)
+Ftip = Init(Ftip,'pressure',3,'friction',0.4);
+deep = false;
+[U,F,SE,eps] = Run(Ftip,deep); % SE = strain energy
+Animate(Ftip,U,F,eps,deep,true)
+N = Ftip.Nl;
ac = contactArea(Ftip,U);
delta = normIndent(Ftip,U);
div = divergence(Ftip,U);
-[sedXY,sedXZ] = strainEnergy(Ftip,U);
+SET = sum(SE([2:N+1,N+3:end],:)); % strain energy total (excluding the bone)
+%[sedXY,sedXZ] = strainEnergy(Ftip,U);
-N = Ftip.Nl;
forceN = abs(sum(F(2:N+1,:),1));
figure;
subplot(1,2,1); hold on
plot(forceN,ac)
ylabel('contact area')
subplot(1,2,2); hold on
-plot(forceN,sedXY)
+plot(forceN,SET)
ylabel('strain energy (mJ/mm3)')
-% save(['Stiffness/' num2str(Ftip.kc/5e-3*1e-3) 'kPa.mat'],'U','F','ac','sed')
-% save(['Friction/' num2str(Ftip.fclb) '.mat'],'U','F','ac','sed')
-% save(['Curvature/' num2str(Ftip.Rcurv) '.mat'],'U','F','ac','sed')
+% save(['../Stiffness/' num2str(Ftip.kc/Ftip.d*2*1e-3) 'kPa.mat'],'U','F','ac','SET')
+% save(['../Friction/' num2str(Ftip.fclb) '.mat'],'U','F','ac','SET')
+% save(['../Curvature/' num2str(Ftip.Rcurv) '.mat'],'U','F','ac','SET')
+