3 views (last 30 days)
Show older comments
K_A on 12 Nov 2021
-
-
Link
Direct link to this question
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary
Commented: Walter Roberson on 28 Nov 2021
Open in MATLAB Online
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
5 Comments Show 3 older commentsHide 3 older comments
Show 3 older commentsHide 3 older comments
Jan on 12 Nov 2021
Direct link to this comment
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary#comment_1830549
I do not understand this sentence: "if particles enter into the cylinder, its previous position should not change".
K_A on 12 Nov 2021
Direct link to this comment
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary#comment_1830579
- cylinder_particle.PNG
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 on 12 Nov 2021
Direct link to this comment
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary#comment_1830594
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 on 12 Nov 2021
Direct link to this comment
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary#comment_1830634
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 on 13 Nov 2021
Direct link to this comment
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary#comment_1832049
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.
Sign in to comment.
Sign in to answer this question.
Answers (1)
Jan on 12 Nov 2021
Edited: Jan on 12 Nov 2021
Open in MATLAB Online
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.
4 Comments Show 2 older commentsHide 2 older comments
Show 2 older commentsHide 2 older comments
K_A on 13 Nov 2021
Direct link to this comment
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary#comment_1831514
Edited: Walter Roberson on 28 Nov 2021
Open in MATLAB Online
- Capture.PNG
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 on 13 Nov 2021
Direct link to this comment
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary#comment_1831904
Open in MATLAB Online
"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 on 13 Nov 2021
Direct link to this comment
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary#comment_1832029
Edited: K_A on 13 Nov 2021
Open in MATLAB Online
- Before_check.PNG
- after_check.png
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 on 28 Nov 2021
Direct link to this comment
https://ch.mathworks.com/matlabcentral/answers/1585289-how-to-stop-particles-crossing-the-solid-cylinder-boundary#comment_1856154
We need your full code, so that we are not busy trying to fix problems that we imagine that you might possibly have.
Sign in to comment.
Sign in to answer this question.
See Also
Categories
MATLABLanguage FundamentalsMatrices and ArraysCreating and Concatenating Matrices
Find more on Creating and Concatenating Matrices in Help Center and File Exchange
Tags
- particle
- indexing
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français
- United Kingdom(English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)
Contact your local office
