function adh = adhesionScheme(o,H, delta, m, n,walls,wallsInt,wallsExt)
% flag whether we have walls or not
confined = 0;
if ~isempty(walls) || ~isempty(wallsInt)
confined = 1;
end
% flag for whether there are two sets of walls with different
% discretization
diffDiscWalls = 0;
if ~isempty(wallsInt)
diffDiscWalls = 1;
end
oc = curve;
% Compute x,y coordinates of net repulsive force on each point of each
% vesicle due to all other vesicles and walls
adh = zeros(2*o.N,o.nv);
[ox,oy] = oc.getXY(o.X);
if confined
if ~diffDiscWalls
[xwalls,ywalls] = oc.getXY(walls.X);
Nbd = walls.N;
else
[xwallsInt,ywallsInt] = oc.getXY(wallsInt.X);
NbdInt = wallsInt.N;
[xwallsExt,ywallsExt] = oc.getXY(wallsExt.X);
NbdExt = wallsExt.N;
end
end
for k = 1:o.nv
adhx = zeros(o.N,1); adhy = zeros(o.N,1);
[~,arcLen] = oc.arcLengthParameter(ox(:,k),oy(:,k));
hmax = max(diff(arcLen));
notk = [(1:k-1) (k+1:o.nv)];
notk_ox = ox(:,notk);
notk_oy = oy(:,notk);
for j = 1 : o.N
if o.nv ~= 1
% Find the distances to each point on a vesicle
dist = ((ox(j,k) - notk_ox).^2 + ...
(oy(j,k) - notk_oy).^2).^0.5;
dF = (delta^(m-n) - dist.^(m-n))./(dist.^(m+2));
adhx(j) = sum(sum((ox(j,k) - notk_ox).*dF));
adhy(j) = sum(sum((oy(j,k) - notk_oy).*dF));
end
if confined
etaWalls = hmax*minDistt;
if ~diffDiscWalls
dist = ((ox(j,k) - xwalls).^2 + (oy(j,k) - ywalls).^2).^0.5;
% Stiffness l^2
dF = (delta.^(m-n) - dist.^(m-n))./dist.^(m+2);
wallAdhx = sum(sum((ox(j,k)-xwalls).*dF)); % why there is multiplication by 2?????
wallAdhy = sum(sum((oy(j,k)-ywalls).*dF)); % Because We calculate adh differently
% Repulsion due to wall has to beat repulsion due to vesicles if a
% vesicle is a sandwich between a vesicle and a wall
adhx(j) = adhx(j) + wallAdhx;
adhy(j) = adhy(j) + wallAdhy;
else
% Interior walls
dist = ((ox(j,k) - xwallsInt).^2 + (oy(j,k) - ywallsInt).^2).^0.5;
L = floor(etaWalls./dist);
% Stiffness l^2
dF = (delta.^(m-n) - dist.^(m-n))./dist.^(m+2);
wallIntAdhx = sum(sum((ox(j,k)-xwallsInt).*dF));
wallIntAdhy = sum(sum((oy(j,k)-ywallsInt).*dF));
% Exterior wall
dist = ((ox(j,k) - xwallsExt).^2 + (oy(j,k) - ywallsExt).^2).^0.5;
L = floor(etaWalls./dist);
% Stiffness l^2
dF = (delta.^(m-n) - dist.^(m-n))./dist.^(m+2);
wallExtAdhx = sum(sum((ox(j,k)-xwallsExt).*dF));
wallExtAdhy = sum(sum((oy(j,k)-ywallsExt).*dF));
% Repulsion due to wall has to beat repulsion due to vesicles if a
% vesicle is a sandwich between a vesicle and a wall
adhx(j) = repx(j) + wallIntAdhx + wallExtAdhx;
adhy(j) = repy(j) + wallIntAdhy + wallExtAdhy;
end % diffDiscWalls
end % confined
end % o.N
% Adhesion on the kth vesicle multiplied with Adhesion Constants
adh(:,k) = -(H * m * delta^n)*[adhx; adhy];
end % o.nv
end % AdhesionScheme
这是以下方程的实现粘附力方程
我们和我的顾问讨论过,看来我的实施没有任何问题。 主要问题是,当我增加点数时,附着力会急剧增加。第二个问题是常数 H 应该是正值,但这会使力变弱并且囊泡不会粘附。当我使用负 H 时,它的效果与文献中的结果几乎完美。
我尝试使用 32 点作为基础来调整值,如下所示:
adj_factor = (o.N / 32);
adh(:,k) = -(H * m * delta^n)*[adhx; adhy]/ adj_factor;
下面只是调试线路,以查看每次运行代码时力如何变化。
disp('Total adhesion force on the first vesicle:');
disp(norm(adh(:,1)));
disp('\n');
total_points = o.N * o.nv;
normalized_adh = adh / total_points;
disp('Normalized adhesion force per point:');
disp(norm(normalized_adh));
当你用黎曼和来近似 γk 上的积分时,你没有忘记 dSy 吗?
然后,增加点数将减小与每个点相关的 dSy 的大小,这应该抵消粘附力的增加。