From d17f10c5cd79ee05d560e8e1194820df3ed21dd4 Mon Sep 17 00:00:00 2001 From: lwillemet <l.willemet@tudelft.nl> Date: Mon, 9 May 2022 08:45:00 +0200 Subject: [PATCH] Viscoelasticity --- Fingertip.m | 193 ++++++++++++++++++++++++++++++++++++++++++++-------- README.md | 14 +++- example.m | 21 +++--- 3 files changed, 190 insertions(+), 38 deletions(-) diff --git a/Fingertip.m b/Fingertip.m index b626037..1eea1c1 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 41e76c2..1950765 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 0e7d317..65514cf 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') + -- GitLab