如何阻止粒子穿过实心圆柱体边界？

Question

我正在做一个随机游走粒子跟踪模拟。在我的问题中，我必须在圆柱体上实现反射边界。在这段代码中，检查意味着如果粒子进入圆柱体，它们之前的位置不应该改变。我的意思是，如果一个粒子越过实体边界，这样的移动就会被拒绝，并且在当前迭代期间，示踪剂会停留在其先前的位置。为澄清起见，我附上了一张图片，您可以在其中看到粒子坐在圆圈中。 但是，颗粒仍然进入气缸。下面是我的代码。

%% loop starts

for i = 1:T            %  T is the total computation time of  simulation
  t = i*dt;               %  dt is the  time step
  %% Particle tracking starts
    
  x_ out = x_previous + u.dt + sqrt(2Ddt);  % x_out--> current position
                                            % x_previous--> previous position
                                            % u -> x-velocity component 
  y_ out = y_previous + v.dt + sqrt(2Ddt);  % y_out--> current position
                                            % y_previous--> previous position
                                            % v -> y-velocity component
    
  %% Check whether particle is inside the cylinder obstacle 
  % x_c and y_c is the center location for cylinder (obstacle) in the flow
  Incircle = sqrt((x_out - x_c).^2 + (y_out - y_c).^2);  
  Particle_enter  = find (Incircle < D/2);    % Particle entered in obstacle 
                                              % having diameter (D)
  %% if particle cross cylinder boundary, do not update its position
  if ismember(Particle_enter,1) 
    kk1 = find(ismember(isIncircle, 1) == 1);
    x_out(kk1) = x_previous(kk1);             % Do not update location
    y_out(kk1) = y_previous(kk1);
  end        
    x_previous = x_out;
    y_previous = y_out;

end

任何修改/更正都会有所帮助 谢谢

Answer 1

由于舍入误差，比较float或double时，使用

a-b > tol  % tol = some error tolerance, say, 1e-10

而不是

a > b

另外，sqrt 会引入一些数值误差，比较平方距离更准确。

尝试将代码中的这些行更改为

Incircle = (x_out - x_c).^2 + (y_out - y_c).^2;  
Particle_enter  = find (abs(Incircle - D^2/4) < tol);

如果你随机生成粒子的初始位置，首先要确保它们不在障碍物内。

最后，如果一个粒子碰到一个反射边界，它不应该被反射，而不是回到原来的位置吗？