How to stop particles crossing the solid cylinder boundary

Hello,

I am doing random walk particle tracking simulation. In my problem, I have to implement reflected boundary on the cylinder. In this code, the check implies that if particles enter into the cylinder, its previous position should not change.

However, particles still entering into the cylinder. Below is my code.

%% 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 of                                                                           
    %  particle, u is the x-velocity component 
    
    y_ out = y_previous + v.dt + sqrt(2Ddt);  % y_out--> current position, y_in--> previous position of                                                                          # 
    %  particle, v is the y-velocity velocity component
    
    %% Check whether particle reflects from cylinder obstacle 
    % x_c and y_c is the center location for cylinder (obstacle) in the flow
    Incircle = sqrt((x_out - x_c(j)).^2 + (y_out - y_c(j)).^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 = x_previous(kk1);                                    % Do not update location
    y_out = y_previous(kk1);
end        
x_previous = x_out;
y_previous = y_out;

end

Any suggestion/modification in the code will be helpful

Thanks

[Jan]

I do not understand this sentence: "if particles enter into the cylinder, its previous position should not change".

[K_A]

Hello Jan,

I mean, if a particle crosses the solid boundary, such a move is rejected, and the tracer stays at its previous position during the current iteration. For clarification, I am have attached an image where you can see particle sitting in the circle.

[Image Analyst]

Try inpolygon(). Wouldn't it be more accurate to have the particle bounce/reflect off the surface rather than stop dead in its tracks?

[K_A]

Hello Image Analyst,

I tried to reflect it from the cylinder boundary by finding those particles which are sitting in the cylinder and kick them

out by giving them reverse velocity. However, that also did not work.

[Image Analyst]

Yes of course -- that would be the wrong formula unless the particle was aiming directly at the center. If it's a "glancing blow" then you need to take the slope of the cylinder at the point of impact into account. You need to have another array, slope, where the elements are the arctangent of y_c over x_c. Then it's just simple geometry/algebra/trigonometry to reflect a line off of a plane at that slope.

Answer 1

[Jan]

What about:

Outcircle = ((x_out - x_c(j)).^2 + (y_out - y_c(j)).^2) >= (D/2)^2;
x_out(Outcircle) = x_previous(Outcircle);
y_out(Outcircle) = y_previous(Outcircle);

Note: sqrt() is expensive. Squaring the radius is much cheaper.

[K_A]

Hello Jan,

I made changes in my code as you suggested. But still it does not work. Particles seems to freeze and they are not moving after applying this check. Result Image is attached (blue is the target cylinder/circle ).

%% 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 of                                                                           
    %  particle, u is the x-velocity component 
    
    y_out = y_previous + v.dt + sqrt(2Ddt);  % y_out--> current position, y_in--> previous position of                                                                          # 
    %  particle, v is the y-velocity velocity component
    Outcircle = ((x_out - x_c(j)).^2 + (y_out - y_c(j)).^2) >= (D/2)^2;
    x_out(Outcircle) = x_previous(Outcircle);
    y_out(Outcircle) = y_previous(Outcircle);
    %% Check whether particle reflects from cylinder obstacle 
    % x_c and y_c is the center location for cylinder (obstacle) in the flow
    plot(x_out, y_out, '.')    
    x_previous = x_out;
    y_previous = y_out;
    
end

[Jan]

"Particles seems to freeze and they are not moving after applying this check."

Yes, of course. This is what you have asked for:

"if particles enter into the cylinder, its previous position should not change."

Please use the tools for format code in the forum. I've done this in your question and in your comment for you. In both cases you see immediately, that there is a strange space:

    x_ out = x_previous + u.dt + sqrt(2Ddt);
    % ^
    y_ out = y_previous + v.dt + sqrt(2Ddt);
    % ^    

This seems to be a typo, but if x_ and y_ are not defined before, it should stop with an error message. So how do you create the output?

In this code you add some fixed values to the x and y positions. Are you sure that you want to access the field dt of the structs u and v? Motion simulations usually have a term like velocity_x * dt to determine the motion of a partical.

Please take into account that "does not work" is not useful to explain a poblem. Do you get an error message? If so, which one? Or do the results differ from your expectations? Then explain both - remember that the readers do not know, what you expect.

[K_A]

Hello Jan,

I apologize for ambiguity and formating of my question. I do store output in x_out and y_out (there is no space in x and out). For clarity, I have modified output variables names to avoid any amibiguties. The physical interpretation of equation

xout = xprevious + u*dt + randn(N,1).*(2*D*dt)^0.5;
yout = yprevious + v*dt + randn(N,1).*(2*D*dt)^0.5;

is that In order to move a particle in the space (here 2-dimension), I need its intial positions (xprevious and yprevious) and corresponding velocities (u--> velocity component in x direction and v--> velocity component in y direction). I obtain velocity components from other solver. In addition, here the walker is random (not purposeful or determinstic), so I add random noise (here N is the total number of particles) with some strength (sqrt(2*D*dt)).

By this equation, particles should move in forward direction. In my problem, there is an obstacle. Particles cannot cross solid boundary, so either they have to reflect or remain where they are. I am trying to implement the later one which is

Outcircle = ((xout - xc).^2 + (yout - yc).^2 >= (D/2)^2);
xout(Outcircle) = xprevious(Outcircle);
yout(Outcircle) = yprevious(Outcircle);

I did not get any error message. When I say it does not work, I mean before applying the check, I got result shown in before_check image. Here particle moves as time advances but enters into the obtscle.

However, when I apply the check as you suggested, I got result shown in after_check image. Here particles are not moving at all which is not expected by the governing equation defined above.

[Walter Roberson]

We need your full code, so that we are not busy trying to fix problems that we imagine that you might possibly have.