GPUSPH / gpusph

The world's first CUDA implementation of Weakly-Compressible Smoothed Particle Hydrodynamics
159 stars 67 forks source link

Rigid body sliding along an inclined plane #12

Open ghost opened 9 years ago

ghost commented 9 years ago

Hi, everyone. I want to simulate a rigid body sliding along an inclined plane. Its motion is influenced by water particles and ramp particles. Can GPUSPH deal with the interaction between the rigid body and the ramp? Thank you!

Narcolessico commented 9 years ago

Hello, GPUSPH relies on the ODE library for body dynamics. You can define a Problem with rigid bodies interacting with planes (see e.g. how dCreatePlane() is used in OdeObject problem). Please note that fluid particles interact with floating bodies but not with ODE planes; to make the fluid interact with the plane you should define the same plane also as GPUSPH plane (see instead fill_planes() and copy_planes() in DamBreak3D problem).

ghost commented 9 years ago

Thank you. And the rigid body could interact with water particles. Is that right? If I have a rigid body sliding along an ramp and then along a flat plane, can I make sure that the rigid body will not penetrate the ramp and the flat plane by making the rigid body interact with the ramp and water?

Narcolessico commented 9 years ago

Yes. Objects (i.e. classes derived from Object) are filled or covered with particles and thus interact automatically with fluid; they also offer methods to declare a correspondent ODE body (ODEBodyCreate() for body dynamics, ODEGeomCreate() for collisions). Planes instead do not have a correspondent class and therefore need separate definitions.

Please note that Problem interface will change in the near future, offering a simpler object definition API and a Plane class. For more clarifications you're welcome to join irc://irc.freenode.net/gpusph

ghost commented 9 years ago

Thanks for your answer. Could the interaction between a plane and a rigid body avoids that the rigid body penetrating the plane? Because I used SPHysics recently and found that it couldn't deal with the interaction between a rigid body and a plane and that the rigid body will penetrate through the plane.

Narcolessico commented 9 years ago

Yes, an ODE plane does not allow other ODE bodies to penetrate. In OdeObjects example, 5 planes are defined to limit the movements of the rigid bodies. In the ODE_near_callback() method you can also handle the number of contact points, bounciness, friction, etc. of the collision.

ghost commented 9 years ago

Thank you very much. So it sounds like GPUSPH is more advanced than SPHysics. Is there any other program that can solve the interaction between a moving rigid body and a plane and at the same time don't allow the rigid body to penetrate the plane? Because I know GPUSPH can only be used on a computer with a Nvidia Graphics card and unfortunately my graphics is intel graphics.

ghost commented 9 years ago

Is there a way that I could ensure that a rigid body only slides along the slope, not rotating?

alda30 commented 9 years ago

As far as I know, high friction ratio (also depends on the geometry of the rigid body, its mass and the slope of the inclined plane) normally makes the rigid body to overturn rather than sliding.

If you need a pure sliding you have to set the friction ratio of the contact properly. I think in ODE you can adjust contact friction by: contact[i].surface.mu = [someValue between 0 to 1]

ghost commented 9 years ago

Thank you. I will try. And I have another question. Someone told me that I could have a pure sliding by setting angular motor. But after I read ODE's user guide, I feel that angular motor is not what I need. Is that right?

alda30 commented 9 years ago

Well, to me it is hard to answer that. It all depends on what you are looking for and how. Some people use physics engines (and SPH softwares) for animation/video Games/ect. purposes, but some other doing engineering simulations. In the first case, only a "visually plausible" behaviour is sought; however; in the latter case you have to set the real (or at least not out-of-range) engineering values for your parameters, and if you set everything well, the observed behaviour would be possibly close to the real one.

To me it seems lowering the friction ratio to a value about 0 would resist the rigid body from overturning. If not, you can even lock the corresponding DOF related to overturning of the rigid body. This seems more easy, and more straightforward rather than using motors.

ghost commented 9 years ago

Thank you. Before you answered my question, what I wanted to know is whether users could lock the corresponding degree of freedom and how to lock it in GPUSPH. But from your answer, I think setting a suitable friction ratio is better because this method is more close to real world.

ghost commented 9 years ago

Closer to the real world*

ghost commented 9 years ago

Hey. Everyone. I have another problem. I read a manual. It tells me that GPUSPH will automatically save results in a data file and an image file. When I use VTKWRITER or TEXTWRITER, there is no image file. And when I use UDPWRITER, there is no image file and data file. Is that manual outdated? Or could I let GPUSPH store data file and image file at the same time?

Azrael3000 commented 9 years ago

I can't comment on the ODE problem as I don't know it very well. The manual that you have read is out of date, so it no longer represents the current code. We are working on a new one but it will take a few months until it is completed. Regarding your specific question, there is no image saved and as far as I'm aware there is no option to save an image file.

ghost commented 9 years ago

Hi, Azrael3000, I am simulating this project. But I meet a few problems.

  1. I create a slope, that is a cube. I rotate it and then use its function method to remove the water under it. However, the result is that the slope is rotated but water is removed before it's rotated. Do you know the reasons
  2. I can compile my case successfully. But when I run, the terminal shows that max. neighbors numbers 137 greater than MAXNEIBSUM 127 and timestep 0 under machine epsilon 1 - requestion quit. Then Simulation end. So basically, my case doesn't run. Do you know possible reasons? Thank you.
Narcolessico commented 9 years ago

Hi,

  1. the rotation is taken into account in Cube::IsInside(), which is used in the unfilling. How do you rotate the cube?
  2. The maximum number of neighbors limits the size of the neighbor list. You can set a higher value in the constructor with m_simparams.maxneibsnum. This can be normal in some topologies or can be the result of a wrong object positioning (e.g. intersecting particles).
ghost commented 9 years ago

Thank you!

  1. I created an Eulerparameter first, called ep. Then I initialize the slope by slope=Cube(Point(r0, r0, r0), lx, ly, lz, ep). After I fill fluid, I use slope.Unfill(parts, r0).
  2. I will first try set m_simparams.maxneibsnum to a larger number to test whether it is too small.
ghost commented 9 years ago

Hey. This time, I set m_simparams.maxneibsnum to a larger number. However, the terminal shows that FATAL: timestep 0 under machine epsilon at iteration 1 - requestion quit..... Is this problem caused by wrong rotation or unfill?

Narcolessico commented 9 years ago

It could be. You could open the simulation output (e.g. with Paraview) and check if there are any overlapping geometries (esp. fluid particles too close to boundary particles).

ghost commented 9 years ago

Well, is there any problem in rotating the slope and slope.Unfill ? The output shows that the slope is rotated but water is removed before the slope is rotated.

ghost commented 9 years ago

Hi, I read Cube::IsInside(), and I find that it contains a variable m_ODERot. I guess this variable should be related to EulerParameter because I use a EulerParameter to rotate a cube. However, EulerParameter only contains a variable called m_rot. And I couldn't find where m_ODERot is defined. So I guess this is the reason that the slope couldn't be rotated. Do you know where the m_ODERot is defined and it's meaning? Thanks.

Narcolessico commented 9 years ago

Do you pass the EulerParameters in the constructor of Cube? Could you show your code (e.g. in a pastebin)?

ghost commented 9 years ago

Wow. So happy. Here is header file.http://pastebin.ca/2975812 Here is .cc file. http://pastebin.ca/2975811

Narcolessico commented 9 years ago

This is a bug in the interface and has been already fixed in a developer branch which will be released hopefully in a few weeks. As a workaround, please try replacing Cube.cc line 423

Point lp = (p - m_origin).TransposeRot(m_ODERot);

with

Point lp = m_ep.TransposeRot(p - m_origin);

then the unfilling should work as expected.

There are minor spacing issues in the problem, e.g. the rigid body touches the shore while it should be r0 far from it. I suggest you to add one object at a time and verify the output before adding anything else (e.g. disable all planes and rigidbody and test first the fluid).

ghost commented 9 years ago

Okay. Thank you sooooo much!

ghost commented 9 years ago

Hi, Narcolessico. I have made some corrections to my case. But new problem appear. When I run my case, terminal always shows that "collision between cube and plane 0 contact points". At first, there are 2 points. Do you know possible reasons? Here are the modified source codes. http://pastebin.ca/2976498

Narcolessico commented 9 years ago

Hello,

did you check the output? What could happen is that the cube penetrates the plane a little already in the initial scenario, and thus it is violently expelled. Try to increase the writer frequency and see how the cube moves in time. Also, note that GPUSPH and ODE use a different equation of plane (ax+by+cz+d=0 vs ax+by+cz=d).

Please use separate threads for separate problems. For generic requests for help, I would suggest you to use the IRC channel instead: irc://irc.freenode.net/gpusph

ghost commented 9 years ago

Yeah. But I could not see planes in Paraview. How do you see planes? And I think GPUSPH and ODE use same equation of plane.

ghost commented 9 years ago

Sorry, I meant that I did think GPUSPH and ODE use same equation of plane before you told me.

Narcolessico commented 9 years ago

I'm afraid you have to manually add it in paraview (sources / plane / insert coords of 3 points).

ghost commented 9 years ago

Could you explain more explicitly? I am not familiar with Paraview.