matlab - 乳液在滚筒之间分流

Question

这是我的问题的真正简化：

我最初有许多相互连接的滚筒，其中一个上带有墨水，并且它们彼此相对旋转。当辊子旋转时，墨水在辊子之间以 50/50 的比例共享。一些辊子有不止一个连接，因此发生分墨的顺序很重要。这是一个可视化：

我基本上是在寻找允许程序按顺序遍历连接的东西。我曾想过这样的事情：

for each roller
    if number_of_connections > 2  % # more than one connection for roller
        % # Iterate through connections sequentially, depending on the direction
        % # of rotation
    end
end

Answer 1

我一直想知道你最终会为这个问题提供什么样的解决方案。从您的评论中：

我现在想的方式是在一个滚轮上（可能是最大的）有一个固定的段分辨率，所有其他的应该是一个相对的、整数的比率，以便段在系统旋转时匹配。

看来您目前正在研究与我在这里拥有的相同类型的解决方案。

注意：最后有一个易于复制的版本。

一开始，我使用了 EitanT 在他对上一个问题的回答中使用的结构。

% # Initial state
C = [0, 0;   % # Roller centers (x, y)
     2, 0;
     2, 4;
     2,-5;
     8, 4;
     8,-5;
     8,-1];
R = [1,1,3,4,3,2,2];    % # Roller radii (r)
N = numel(R);           % # Amount of rollers

% # Draw the rollers
figure, hold on
ang = 0:0.1:(2 * pi);
for i = 1:N
    plot(C(i, 1) + R(i) * cos(ang), C(i, 2) + R(i) * sin(ang))
    text(C(i, 1), C(i, 2), num2str(i))
end
title('Ink rollers'), axis image

% # Find connected rollers
isconn = @(m, n)(sum(([1, -1] * C([m, n], :)) .^ 2) - sum(R([m, n])) .^ 2 < eps);
[Y, X] = meshgrid(1:N, 1:N);
conn = reshape(arrayfun(isconn, X(:), Y(:)), N, N) - eye(N);

我使用此示例查看连接顺序的差异。这是生成的滚筒图像：

这个想法是第一辊充满墨水，并且这种墨水通过系统传播通过辊 3 和 5 或 4 和 6 到达辊 7。

我手动设置了所有滚轮的方向。这一步可以通过使用 connection-matrix 自动完成conn。

% # Direction of rotation (clockwise = -1, anticlockwise = 1)
rotDir = [1,-1,1,1,-1,-1,1];

然后，我指定最小滚筒的箱数，并按比例扩大所有其他滚筒的箱数。

% # Number of bins specified for smallest roller
nBins_min = 20;
nBins = round(nBins_min*R/min(R));

下一步是初始化滚筒。我使用了一个结构将墨水、连接和滚轮方向保存在同一个变量中。这个想法是保持墨水值并跟踪所有辊子的每个部分的连接。如果jj滚子ii的一部分未连接到另一个滚子，则用 中的零表示rollers(ii).connections(jj)。否则，如果它已连接，则此单元元素将包含它所连接的滚轮的滚轮索引。

% # Initialize roller struct
rollers = struct('ink',{},'connections',{},'rotDirection',{});

% # Ink
for ii = 1:N
    rollers(ii).ink = zeros(1,nBins(ii));
end
rollers(1).ink = ones(1,nBins(1));

% # Connections
for ii = 1:N
    rollers(ii).connections = zeros(1,nBins(ii));
end
for ii = 1:N
    for jj = 1:N
        if(ii~=jj)
            if(conn(ii,jj) == 1)
                connInd = getConnectionIndex(C,ii,jj,nBins(ii));
                rollers(ii).connections(connInd) = jj;
            end
        end
    end
end

% # Direction of rotation
for ii = 1:N
    rollers(ii).rotDirection = rotDir(ii);
end

我getConnectionIndex()以以下（非常丑陋）的方式实现了上面使用的功能：

function connectionIndex = getConnectionIndex(C,ii,jj,nBins)

p1 = C(ii, :);
p2 = C(jj, :);

if(abs(p1(2)-p2(2))<eps)
    if(p2(1)>p1(1))
        angle = 0;
    else
        angle = pi;
    end
elseif(abs(p1(1)-p2(1))<eps)
    if(p2(2)>p1(2))
        angle = pi/2;
    else
        angle = 3*pi/2;
    end
else
    angle = mod( atan((p2(1)-p1(1))/(p2(2)-p1(2))), 2*pi);
end

connectionIndex = 1+floor(nBins*angle/(2*pi));

end

该函数使用滚轮的中心点来获取连接的相应角度。从这个角度值计算段的索引。

每个时间步长（所有段旋转一个步长）计算并保存每个辊上的平均墨水量。初始化这个矩阵，计算并保存初始墨水分布。

% # Initialize averageAmountOfInk and calculate initial distribution
nTimeSteps = 200;
averageAmountOfInk = zeros(nTimeSteps,N);
for ii = 1:N
    averageAmountOfInk(1,ii) = mean(rollers(ii).ink);
end

对每个时间步执行以下步骤。

% # Iterate through timesteps
for tt = 2:nTimeSteps

第一个辊子充满墨水，所有辊子根据其旋转方向旋转一步。

% # Fill first roller with ink
rollers(1).ink = ones(1,nBins(1));

% # Rotate all rollers
for ii = 1:N
        rollers(ii).ink(:) = ...
            circshift(rollers(ii).ink(:),rollers(ii).rotDirection);
end

然后通过找到匹配的连接并将这两个连接的墨水平均分配来更新所有辊连接。

% # Update all roller-connections
for ii = 1:N
    for jj = 1:nBins(ii)
        if(rollers(ii).connections(jj) ~= 0)
            index1 = rollers(ii).connections(jj);
            index2 = find(ii == rollers(index1).connections);
            ink1 = rollers(ii).ink(jj);
            ink2 = rollers(index1).ink(index2);
            rollers(ii).ink(jj) = (ink1+ink2)/2;
            rollers(index1).ink(index2) = (ink1+ink2)/2;
        end
    end
end

作为最后一步，计算辊上的平均墨水量，并在循环完成后绘制这些值。

    % # Calculate average amount of ink on each roller
    for ii = 1:N
        averageAmountOfInk(tt,ii) = mean(rollers(ii).ink);
    end
end

figure
plot(averageAmountOfInk,'b')
xlabel('Timesteps')
ylabel('Ink')

运行代码，它会生成以下所有滚筒的平均墨水图。

对于最小滚筒的 20 段和 60 个时间步，我们得到下图：

如果我们运行 2000 个时间步长的模拟，我们可以看到收敛到一个，因为所有的滚筒都充满了墨水。

易于复制的版本：

function averageAmountOfInk = inkRollerModel()
% # Initial state
C = [0, 0;   % # Roller centers (x, y)
     2, 0;
     2, 4;
     2,-5;
     8, 4;
     8,-5;
     8,-1];
R = [1,1,3,4,3,2,2];    % # Roller radii (r)
N = numel(R);           % # Amount of rollers

% # Draw the rollers
figure, hold on
ang = 0:0.1:(2 * pi);
for i = 1:N
    plot(C(i, 1) + R(i) * cos(ang), C(i, 2) + R(i) * sin(ang))
    text(C(i, 1), C(i, 2), num2str(i))
end
title('Ink rollers'), axis image

% # Find connected rollers
isconn = @(m, n)(sum(([1, -1] * C([m, n], :)) .^ 2) - sum(R([m, n])) .^ 2 < eps);
[Y, X] = meshgrid(1:N, 1:N);
conn = reshape(arrayfun(isconn, X(:), Y(:)), N, N) - eye(N);

% # Direction of rotation (clockwise = -1, anticlockwise = 1)
rotDir = [1,-1,1,1,-1,-1,1];

% # Number of bins for smallest roller
nBins_min = 20;
nBins = round(nBins_min*R/min(R));

% # Initialize roller struct
rollers = struct('ink',{},'connections',{},'rotDirection',{});

% # Ink
for ii = 1:N
    rollers(ii).ink = zeros(1,nBins(ii));
end
rollers(1).ink = ones(1,nBins(1));

% # Connections
for ii = 1:N
    rollers(ii).connections = zeros(1,nBins(ii));
end
for ii = 1:N
    for jj = 1:N
        if(ii~=jj)
            if(conn(ii,jj) == 1)
                connInd = getConnectionIndex(C,ii,jj,nBins(ii));
                rollers(ii).connections(connInd) = jj;
            end
        end
    end
end

% # Direction of rotation
for ii = 1:N
    rollers(ii).rotDirection = rotDir(ii);
end

% # Initialize averageAmountOfInk and calculate initial distribution
nTimeSteps = 200;
averageAmountOfInk = zeros(nTimeSteps,N);
for ii = 1:N
    averageAmountOfInk(1,ii) = mean(rollers(ii).ink);
end

% # Iterate through timesteps
for tt = 2:nTimeSteps
    % # Fill first roller with ink
    rollers(1).ink = ones(1,nBins(1));

    % # Rotate all rollers
    for ii = 1:N
            rollers(ii).ink(:) = ...
                circshift(rollers(ii).ink(:),rollers(ii).rotDirection);
    end

    % # Update all roller-connections
    for ii = 1:N
        for jj = 1:nBins(ii)
            if(rollers(ii).connections(jj) ~= 0)
                index1 = rollers(ii).connections(jj);
                index2 = find(ii == rollers(index1).connections);
                ink1 = rollers(ii).ink(jj);
                ink2 = rollers(index1).ink(index2);
                rollers(ii).ink(jj) = (ink1+ink2)/2;
                rollers(index1).ink(index2) = (ink1+ink2)/2;
            end
        end
    end

    % # Calculate average amount of ink on each roller
    for ii = 1:N
        averageAmountOfInk(tt,ii) = mean(rollers(ii).ink);
    end
end

figure
plot(averageAmountOfInk,'b')
xlabel('Timesteps')
ylabel('Ink')

end

function connectionIndex = getConnectionIndex(C,ii,jj,nBins)

p1 = C(ii, :);
p2 = C(jj, :);

if(abs(p1(2)-p2(2))<eps)
    if(p2(1)>p1(1))
        angle = 0;
    else
        angle = pi;
    end
elseif(abs(p1(1)-p2(1))<eps)
    if(p2(2)>p1(2))
        angle = pi/2;
    else
        angle = 3*pi/2;
    end
else
    angle = mod( atan((p2(1)-p1(1))/(p2(2)-p1(2))), 2*pi);
end

connectionIndex = 1+floor(nBins*angle/(2*pi));

end

编辑：原始问题值

C = [-276.4, 565.08;...  % # Duct
     -27.82, 616.11;...  % # r2
      41.26, 562.41;...  % # r3
      52.12, 473.07;...  % # r4
     -44.97, 366.25;...  % # ink drum
      137.22, 443.76;... % # r6
      99.32, 362.13;...  % # r7
      141.22, 272.79;... % # r8
      51.67, 237.7;...   % # r9
      173.99, 177.07;... % # r10
     -203.02, 230.52;... % # r11
     -110.9, 213.53;...  % # r12
     -207.33, 131.94;... % # r13
     -187.4, 330.49;...  % # r14
      0,0...             % # Plate cylinder
];% # Roller centres (x, y)

R = [...
95/2,...  % # Duct
80/2,...  % # r2
95/2,...  % # r3
85/2,...  % # r4
208/2,... % # Drum
96/2,...  % # r6
85/2,...  % # r7
112.35/2,... % # r8
81/2,...  % # r9
90/2,...  % # r10
112.35/2,... % # r11
75/2,...  % # r12
86/2,...  % # r13
90/2,...  % # r14
406.5/2 % # Plate
]; % # Roller radii (r)

Answer 2

这实际上是评论，而不是答案。

您和 EitanT 似乎正在融合的方法似乎是一个有趣的数学问题，可以转化为计算机代码。但我有点怀疑你的数学模型是否会导致对真实物理系统行为的准确预测。我能想到的一些问题：

1. 您谈到了应用墨水拆分的顺序。这种方式对于首次通过系统是有意义的，但是一旦所有辊子都有墨水，我就不清楚是否有任何合理的“顺序”可以迭代辊子。如果最初有多个辊子有墨水，那么即使是第一次通过也确实没有意义。
2. 较小的滚子比较大的滚子转动得更快，因此使算法的单次迭代对应于所有滚子的单次旋转实际上并不奏效。
3. 墨水在辊子之间不断地散布，因此每个辊子上的墨水高度是辊子上位置的函数。

换句话说，你似乎忽略了物理系统中的许多复杂性。也许这没问题，这种方法可以捕捉 90% 的重要内容。这会让我感到惊讶……但有趣的是，如此复杂的系统可以如此轻松地建模！

添加：

我刚刚意识到，在您关于此主题的其他问题中，您循环通过滚筒的顺序并不重要。所以“分裂的顺序”在这里可能是一个红鲱鱼。

我也确信您可以从简单的矩阵乘法中得到结果，但还没有尝试计算出相关矩阵。但这可能很容易弄清楚！