Commits

Mike Graham committed e5ce04e Draft

Add the old Matlab code.

  • Participants
  • Parent commits 6a9ada8

Comments (0)

Files changed (4)

matlab/AttardElasticModel.m

+%% ELASTIC ATTARD CONTACT MODEL
+%   This script calculates the interaction between to bodies in contact
+%   based on realistic surface forces. The approach is based on the
+%   work of Attard and Parker (1992) and Attard (2000, 2001).
+%
+%  Author:
+%    Mike Graham
+%    mikegraham@gmail.com
+
+
+%%
+clear
+clc
+
+
+%% SET THE VALUE OF PARAMETERS
+
+% MATERIAL PARAMETERS
+%
+%  Effective modulus
+E = 100 * 10^(6-18); %N/nm^2  =  Exapascals = 10^18 Pa
+%
+%
+% Here is the pressure function p(h).
+%==Exponential-Expoential Potential==
+% P_attract = 12 * 10^(6-18); %Exapascals
+% kappa_attract = 1.0; % 1/(nm)
+% P_repulse = 60 * 10^(6-18); %Exapascals
+% kappa_repulse = 2.0; % 1/(nm)
+%
+% pressure = @(h) P_repulse * exp(-kappa_repulse*h) ...
+%               - P_attract * exp(-kappa_attract*h);
+%
+%==Lennard-Jones adhesion==
+A = 10^(-19 + 9); %nJ = N * nm, Hammicker's constant
+z0 = 1.0; %nm
+pressure = @(h) A/6/pi./h.^3 .* ((z0 ./ h).^6 - 1);
+%
+%==Exponential repulsion==
+% Pr = 50 * 10^(6-18); %EPa
+% kr = 1; %nm
+% pressure = @(h)  Pr * exp(-kr * h);
+%
+%
+% GEOMETRIC AND SIMULATION PARAMETERS
+%
+% The history of the _position_ h_0 through the simulation at step.
+step = 0.05;
+h0 = [5:-step:-5  -5:step:5];
+%
+%
+% Specify the profile function
+profilefun = @(r) roundcone(r, 60, 1);
+%
+% The maximum radius to consider (in place of infinity. The point at
+% which no deformation is expected).
+maximum_radius = 25; %nm
+%
+% The number of grid points for the spatial integration.
+meshnum = 400;
+%
+% The maxiumum blending parameter for self-consistent blending. It is
+% the greatest proportion to which the new trial solution is introduced.
+% Smaller is more stable.
+maxblend = 0.01;
+%
+% The error tolerance, which is the largest RMS value that two trials'
+% predictions for the new deformation profile can vary.
+maxerror = 0.001;  %1pm
+
+
+%% INITIALISE AND CALCULATE NEEEDED VARIABLES
+%
+% Set up the time parameters.
+num_time_steps = length(h0);
+%
+%
+% This just allocates a lot of arrays, for speed.
+u = zeros(meshnum, num_time_steps);
+u_inf = zeros(meshnum, num_time_steps);
+pdot = zeros(meshnum, num_time_steps);
+udot = zeros(meshnum, num_time_steps);
+h = zeros(meshnum, num_time_steps);
+integralp = zeros(meshnum, 1);
+integralpdot = zeros(meshnum, 1);
+udotnew = zeros(meshnum, 1);
+F = zeros(1, num_time_steps);
+u0 = zeros(1, num_time_steps);
+%
+% Make the points in the radii.
+radius = linspace(0, maximum_radius, meshnum);
+dr = radius(2);
+radius = radius + dr/2;
+%
+% This is the profile for the undeformed shape, i.e. round (for radius
+% values much smaller than the effective radius of the solids.)
+profile = profilefun(radius');
+%
+% The initial separation should be the initial undeformed separation.
+h(:,1) = h0(1) + profile;
+%
+Deltah0 = diff(h0);
+%
+Delta_u0_trial = 0;
+
+
+%% COMPUTE THE KERNEL
+% The kernel is used to perform the pressure-force relation (Attard 2001,
+% Eq. 2) and arises from the enforcement of axisymmetry
+% (Attard and Parker 1992, 15).
+kernel = zeros(meshnum, meshnum);
+for r = 1:meshnum
+for s = 1:meshnum
+  % For the special case that the dummy variable equals the radius value,
+  % the elliptic integral is unbounded and we have to calculate the
+  % kernel using the limit from each side (Attard 2001).
+  if s==r
+    Iminus = (dr/pi - dr^2/4/pi/radius(r)) ...
+             * (1 + log(16*radius(r)^2/(radius(r)*dr - dr^2/4)) );
+    Iplus  = dr/4*log(16*(radius(r)+dr/2)^2/(radius(r)*dr + dr^2/4))...
+             + 4*radius(r)/pi*log((2*radius(r)+dr)/(2*radius(r)+dr/2));
+    kernel(r,s) = (Iminus+Iplus)/radius(r)/dr;
+    break
+  else
+    % Compute the kernel normally. Since it is symmetric k(r,s) = k(s,r),
+    % we only do half of it and the break above keeps us from seeing the
+    % case s>r.
+    kernel(r,s) = 4/pi/radius(r) * ellipke((radius(s)/radius(r))^2);
+    kernel(s,r) = kernel(r,s);
+  end
+end
+end
+
+
+%% MAIN SIMULATION
+for t = 2:num_time_steps
+  err = Inf;
+
+  u(:,t) = 0.5 * ( (h0(t) + profile) - h(:,t-1) );
+
+  c = 0;
+
+  blend = maxblend;
+
+  while err > maxerror
+
+    % Calculate the separation (for finding the pressure and pressure
+    % derivative), based on Attard and Parker (1992) Eqs. (17) and (18).
+    h(:,t) = h0(t) + profile - 2 * u(:,t);
+
+    % This is the main equation! (Attard 2000)
+    %   u(r) = int( p(s) k(r,s) s ds, 0, inf)
+    for r = 1:meshnum
+      integralp(r) = sum(pressure(h(:,t)) .* kernel(:,r) .* radius(:)) * dr;
+    end
+    unew  = -2/E * integralp(:);
+
+    % Self-consistent blending:
+    u_old = u(:,t);
+    u(:,t) = (1 - blend) * u(:,t) + blend * unew(:);
+    u0(t) = (3 * u(1,t) - u(2,t)) / 2;
+
+    % Calculate the error
+    if c > 5000
+      error('Convergence judged unlikely')
+    elseif c > 0
+      err = norm(u(:,t) - unew) / meshnum;
+    end
+    [t h0(t) err c/100]
+    c = c + 1;
+
+    % Decrease the blending ratio if we have oscillation or error
+    % growth (Attard 2000).
+    Delta_u0_old = Delta_u0_trial;
+    Delta_u0_trial = u0(t) - u0(t-1);
+    if (Delta_u0_trial * Delta_u0_old) < 0
+      blend = blend / 1.5;
+%    elseif abs(Delta_u0_trial) > abs(Delta_u0_old)
+%      blend = blend / 1.2;
+    elseif c > 0
+      blend = min(1.3 * blend, maxblend);
+    end
+
+%     subplot(2,2,4,'replace');
+%       plot(radius, h(:,t), 'r', radius,h0(t)+profile, 'b:')
+%       title('Rigid and Deformed profiles')
+%       xlabel('Radial Coordinate')
+%       ylabel('Separation, Position')
+%       axis([0 500 0 12])
+
+  end
+
+  % Integrate the pressure to find the load. [Attard and Parker (1992), (20)].
+  %   F = 2pi int(p(r) r dr, r=0, r=inf)
+  F(t) = 2 * pi * sum(pressure(h(:,t)) .* radius(:)) * dr;
+
+  figure(1)
+  subplot(2,2,1,'replace')
+    plot(radius, h(:,t), 'r', radius,h0(t)+profile, 'b:')
+    title('Rigid and Deformed profiles')
+    xlabel('Radial Coordinate')
+    ylabel('Separation, Position')
+    axis([0 10 0 12])
+  subplot(2,2,2,'replace');
+    plot(radius,h(:,1:t));
+    title('h');
+    set(get(gca,'YLabel'),'String','separation (nm)');
+    set(get(gca,'XLabel'),'String','radial coordinate (nm)');
+  subplot(2,2,3,'replace');
+    plot(radius,u(:,1:t));
+    title('u');
+    set(get(gca,'YLabel'),'String','displacement (nm)');
+    set(get(gca,'XLabel'),'String','radial coordinate (nm)');
+  subplot(2,2,4,'replace');
+    plot(radius,pressure(h(:,1:t)));
+    title('p');
+    set(get(gca,'YLabel'),'String','pressure (MPa)');
+    set(get(gca,'XLabel'),'String','radial coordinate (nm)');
+
+  figure(2)
+  subplot(2,2,1,'replace');
+    plot(F(1:t),-u(1,1:t));
+    title('F-\delta');
+    set(get(gca,'XLabel'),'String','force (N)');
+    set(get(gca,'YLabel'),'String','center displacement (nm)');
+  subplot(2,2,2,'replace');
+    plot(h0(1:t),-u(1,1:t));
+    title('deflection vs. position');
+    set(get(gca,'XLabel'),'String','position (nm)');
+    set(get(gca,'YLabel'),'String','center displacement (nm)');
+  subplot(2,2,3,'replace');
+    plot(h0(1:t),F(1:t));
+    title('Force vs. position');set(get(gca,'XLabel'),'String','position (nm)')
+    set(get(gca,'YLabel'),'String','force (N)');
+end
+delta = -2*u(1,:);
+
+
+%% REFERENCES
+% Attard, P., and Parker, J. L. (1992). "Deformation and adhesion of elastic
+%   bodies in contact." _Physical Review A_, 46(12), 7959-7971
+%
+% Attard, P. (2000). "The Interaction and Deformation of Elastic Bodies The
+%   Origin of Adhesion Hysteresis". _Journal of Physical Chemistry_, 104,
+%   10635-10641
+%
+% Attard, P. (2001). "Interaction and deformation of viscoelastic particles:
+%   Nonadhesive particles." _Physical Review E_, 63, 1-9

matlab/AttardViscoelasticModel.m

+%% VISCOELASTIC ATTARD CONTACT MODEL
+%   This script calculates the interaction between two viscoelastic bodies in
+%   contact based on realistic surface forces. The approach is based on the
+%   work of Attard and Parker (1992) and Attard (2000, 2001a, b).
+%
+%  Author:
+%    Mike Graham
+%    mikegraham@gmail.com
+
+%%
+clear
+clc
+
+
+%% SET THE VALUE OF PARAMETERS
+
+% MATERIAL PARAMETERS
+%
+%  Effective modulus
+E_0 = 200 * 10^(9-18); %N/nm^2  =  Exapascals = 10^18 Pa
+E_inf = 100 * 10^(9-18); %N/nm^2  =  Exapascals = 10^18 Pa
+tau = 0.5;
+%
+%
+% Here is the pressure function p(h).
+%==Exponential-Expoential Potential==
+% P_attract = 12 * 10^(6-18); %Exapascals
+% kappa_attract = 1.0; % 1/(nm)
+% P_repulse = 60 * 10^(6-18); %Exapascals
+% kappa_repulse = 2.0; % 1/(nm)
+%
+% pressure = @(h) P_repulse * exp(-kappa_repulse*h) ...
+%               - P_attract * exp(-kappa_attract*h);
+%
+%==Lennard-Jones adhesion==
+% A = 10^(-19 + 9); %nJ = N * nm, Hammicker's constant
+% z0 = 1.0; %nm
+% pressure = @(h) A/6/pi./h.^3 .* ((z0 ./ h).^6 - 1);
+%
+%==Poisson-Boltzman repulsion==
+Pr = 100 * 10^(12-18); %EPa
+kr = 3.0; %1/nm
+pressure = @(h)  Pr * exp(-kr * h);
+
+% GEOMETRIC AND SIMULATION PARAMETERS
+%
+% The history of the "rigid seperation" h_0 (at the center) through the
+% simulation at each time.
+step = 0.01;
+h0 = [10:-step:-5  -5:step:10];
+time = linspace(0, 1, length(h0));
+%
+%
+% Specify the profile function
+profile = @(r) roundcone(r, 60, 2);
+%
+% The maximum radius to consider (practically infinity, at which no deformation
+% is expected).
+maximum_radius = 25; %nm
+%
+% The number of grid points for the spatial integration.
+meshnum = 600;
+%
+% The blending parameter for self-consistent blending. It is the proportion to
+% which the new trial solution is introduced.
+maxblend = 0.15;
+%
+% The error tolerance, which is the largest RMS value that two trials'
+% predictions for the new deformation profile can vary.
+maxerror = 0.0005;
+
+
+%% INITIALISE AND CALCULATE NEEEDED VARIABLES
+
+% Set up the time parameters.
+num_time_steps = length(time);
+time_step = diff(time);
+
+% This just allocates a lot of arrays, for speed.
+u = zeros(meshnum, num_time_steps);
+u_inf = zeros(meshnum, num_time_steps);
+pdot = zeros(meshnum, num_time_steps);
+integralp = zeros(meshnum, 1);
+integralpdot = zeros(meshnum, 1);
+udotnew = zeros(meshnum, 1);
+udotold = zeros(meshnum, 1);
+udot = zeros(meshnum, num_time_steps);
+h = zeros(meshnum, num_time_steps);
+F = zeros(1, num_time_steps);
+u0 = zeros(1, num_time_steps);
+
+% Make the points in the radii.
+radius = linspace(0, maximum_radius, meshnum);
+dr = radius(2);
+radius = radius + dr/2;
+
+% This is the profile for the undeformed shape, i.e. round (for radius values
+% much smaller than the effective radius of the solids.)
+profile = profile(radius');
+
+% The initial separation should be the initial undeformed separation.
+h(:,1) = h0(1) + profile;
+
+Deltah0 = diff(h0);
+
+
+
+Delta_u0_trial = 0;
+
+
+%% COMPUTE THE KERNEL
+% The kernel is used to perform the pressure-force relation (Attard 2001a, 2)
+% and arises from the enforcement of axisymmetry (Attard and Parker
+% 1992, 15).
+kernel = zeros(meshnum, meshnum);
+for r = 1:meshnum
+for s = 1:meshnum
+  % For the special case that the dummy variable equals the radius value, the
+  % elliptic integral is unbounded and we have to calculate the kernel using the
+  % limit from each side (Attard 2001a).
+  if s==r
+    Iminus = (dr/pi - dr^2/4/pi/radius(r)) ...
+             * (1 + log(16*radius(r)^2/(radius(r)*dr - dr^2/4)) );
+    Iplus  = dr/4*log(16*(radius(r)+dr/2)^2/(radius(r)*dr + dr^2/4))...
+             + 4*radius(r)/pi*log((2*radius(r)+dr)/(2*radius(r)+dr/2));
+    kernel(r,s) = (Iminus+Iplus)/radius(r)/dr;
+    break
+  else
+    kernel(r,s) = 4/pi/radius(r) * ellipke((radius(s)/radius(r))^2);
+    kernel(s,r) = kernel(r,s);
+  end
+end
+end
+
+
+%% MAIN SIMULATION
+for t = 2:num_time_steps
+  err = Inf;
+
+  u(:,t) = h0(t) + profile - h(:,t-1);
+
+  c = 0;
+
+  blend = maxblend;
+
+  while err > maxerror
+    % Calculate the separation (for finding the pressure and pressure
+    % derivative), based on Attard and Parker (1992) equations (17) and (18).
+    h(:,t) = h0(t) + profile - u(:,t);
+
+    for r = 1:meshnum
+      pressuredot(r) = (pressure(h(r,t)) - pressure(h(r,t-1))) / time_step(t-1);
+      %pressuredot(r) = pressureprime(h(r,t)) * h0-udot
+    end
+
+    for r = 1:meshnum
+      integralpdot(r) = sum(pressuredot(:) .* kernel(:,r) .* radius(:)) * dr;
+    end
+
+    for r = 1:meshnum
+      integralp(r) = sum(pressure(h(:,t)) .* kernel(:,r) .* radius(:)) * dr;
+    end
+    u_inf  = -2/E_inf * integralp(:);
+
+
+    % Attard (2001a) (5)
+    for r = 1:meshnum
+      udot(r,t) = -(1/tau)*(u(r,t) - u_inf(r)) - 2/E_0*integralpdot(r);
+    end
+
+
+    udotnew = (1 - blend) * udotold(:) + blend * udot(:,t);
+
+    %udotold = udot(:,t);
+
+    % Self-consistent blending:
+
+%    u(:,t) = (1 - blend) * u(:,t) + blend * unew(:);
+%    u0(t) = (3 * u(1,t) - u(2,t)) / 2;
+
+    % Calculate the error
+    if c > 10000
+      error('Convergence judged unlikely')
+    elseif c > 30
+      err = norm(udotnew - udotold) *  time_step(t-1) / meshnum;
+    end
+    [t h0(t) err c/100]
+    c = c + 1;
+
+    udotold = udotnew;
+    u(:,t) = u(:,t-1) + udotnew * time_step(t-1);
+
+
+    % Decrease the blending ratio if we have oscillation or error growth (Attard 2000).
+    Delta_u0_old = Delta_u0_trial;
+    Delta_u0_trial = u0(t) - u0(t-1);
+    if (Delta_u0_trial * Delta_u0_old) < 0
+      blend = blend / 1.5;
+%    elseif abs(Delta_u0_trial) > abs(Delta_u0_old)
+%      blend = blend / 1.2;
+    elseif c > 0
+      blend = min(1.3 * blend, maxblend);
+    end
+
+%     subplot(2,2,4,'replace');
+%     plot(radius, h(:,t), 'r', radius,h0(t)+profile, 'b:')
+%     title('Rigid and Deformed profiles')
+%     xlabel('Radial Coordinate')
+%     ylabel('Separation, Position')
+%     axis([0 500 0 12])
+
+  end
+
+
+  % Integrate the pressure to find the load. [Attard and Parker (1992), (20).]
+  %   F = 2pi int(p(r) r dr, 0, inf)
+  F(t) = 2 * pi * sum(pressure(h(:,t)) .* radius(:)) * dr;
+
+
+%   figure(1)
+%   subplot(2,2,1,'replace')
+%     plot(radius, h(:,t), 'r', radius,h0(t)+profile, 'b:')
+%     title('Rigid and Deformed profiles')
+%     xlabel('Radial Coordinate')
+%     ylabel('Separation, Position')
+%     axis([0 15 0 20])
+%   subplot(2,2,2,'replace');
+%     plot(radius,h(:,1:t));
+%     title('h');
+%     set(get(gca,'YLabel'),'String','separation (nm)');
+%     set(get(gca,'XLabel'),'String','radial coordinate (nm)');
+%   subplot(2,2,3,'replace');
+%     plot(radius,u(:,1:t));
+%     title('u');
+%     set(get(gca,'YLabel'),'String','displacement (nm)');
+%     set(get(gca,'XLabel'),'String','radial coordinate (nm)');
+%   subplot(2,2,4,'replace');
+%     plot(radius,pressure(h(:,1:t)));
+%     title('p');
+%     set(get(gca,'YLabel'),'String','pressure (MPa)');
+%     set(get(gca,'XLabel'),'String','radial coordinate (nm)');
+%
+%   figure(2)
+%   subplot(2,2,1,'replace');
+%     plot(F(1:t),-u(1,1:t));
+%     title('F-\delta');
+%     set(get(gca,'XLabel'),'String','force (N)');
+%     set(get(gca,'YLabel'),'String','center displacement (nm)');
+%   subplot(2,2,2,'replace');
+%     plot(h0(1:t),-u(1,1:t));
+%     title('deflection vs. position');
+%     set(get(gca,'XLabel'),'String','position (nm)');
+%     set(get(gca,'YLabel'),'String','center displacement (nm)');
+%   subplot(2,2,3,'replace');
+%     plot(h0(1:t),F(1:t));
+%     title('Force vs. position');set(get(gca,'XLabel'),'String','position (nm)')
+%     set(get(gca,'YLabel'),'String','force (N)');
+end
+delta = -2*u(1,:);
+
+
+tau
+
+%% REFERENCES
+% Attard, P., and Parker, J. L. (1992). "Deformation and adhesion of elastic
+%   bodies in contact." _Physical Review A_, 46(12), 7959-7971
+%
+% Attard, P. (2000). "The Interaction and Deformation of Elastic Bodies The
+%   Origin of Adhesion Hysteresis". _Journal of Physical Chemistry_, 104,
+%   10635-10641
+%
+% Attard, P. (2001a).
+%
+% Attard, P. (2001b).
+These are the old Matlab programs that preceded the 
+Python implementation.

matlab/roundcone.m

+% ROUNDCONE The height of a cone with a rounded tip.
+%    Y = ROUNDCONE(R, ANGLE, RADIUS) returns the value Y of the
+%    height of a body for the a point R away from the centre of
+%    axisymmetry, with an angle ANGLE in degrees and a tip radius
+%    RADIUS.
+%
+%    ROUNDCONE works for R > 0 (or, rather, R > -RADIUS).
+
+% Author:
+%   Mike Graham
+%   mikegraham@gmail.com
+
+function y = roundcone(r, angle, radius)
+  slope = cot(angle/2 * pi/180)
+  r_meet = slope * sqrt(radius^2 / (slope^2 + 1))
+  y_meet = slope * r_meet;
+
+  y_noshift = mycircle(r_meet, radius, 0);
+  shift = y_meet - y_noshift;
+
+  y = slope * r;
+  isround = find(r < r_meet);
+  y(isround) = mycircle(r(isround), radius, shift);
+
+function y = mycircle(r, radius, shift)
+  y = radius - sqrt(radius^2 - r.^2) + shift;