Viscoelasticity

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
         

README.md

 # 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.
 
 ![model_connectivity](https://user-images.githubusercontent.com/58175648/115708505-61f45800-a370-11eb-9294-5a7a5baaf0d0.png)
 
-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.

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')
+