Friction Model

[bencannell]

Hi @tkoziara,

We would like to include functionality to model friction contact within parmec.

What we are trying to achieve is to model some lateral shear resistance between points on particles, which is proportional to the load experienced by a spring which connects these particles. For example, we have springs which model direct vertical contact between vertically stacked components, we would like to model planar friction behaviour between these components. The normal to this plane would be defined by the associated direct contact spring.

Do you think it would be possible to add this into the SPRING definition? Or perhaps a separate command?

Perhaps @stevebrasier, could add some extra clarity to this request?

Best regards,

Ben

[tkoziara]

Most naturally, an optional parameter friction could be added to the SPRING command to accompany the contact plane spring model or a two-point-spring with constant direction. A ratio ks/kn could be used to inherit the tangential spring stiffness from the normal one; similarly, the normal damper curve could be reused in the tangential direction (same ratio could be applied); the normal direction would be the constant direction or the plane normal; by default friction = 0 and defining friction > 0 would enable this; history of local tangential displacements could be tracked and thus stick-slip behaviour could be modelled; Example interface:

SPRING(... , friction, kskn)

Would that be ok?

[bencannell]

Hi @tkoziara,

Yes I think your suggested interface would be appropriate.

Many thanks,

Ben

[tkoziara]

This is now initially implemented and tested; Please see this blog post; closing issue for now.

[bencannell]

Thank you @tkoziara, that was quick! I'll test this week

[bencannell]

Hi @tkoziara,

I've tested the friction definition and I think I may have mislead you regarding the ks/kn definition. I would not expect there to be any tangential stiffness/damping inherited from the normal spring definition, rather just a single tangential force which is the normal spring force * friction coefficient. This force would be applied to the component of the particle's velocity tangential to the normal force (i.e. tangential to the spring direction), and would "effectively" become zero once the velocity of the particle becomes zero.

I've attached a test case which models a particle dropping onto a plane under gravitational loading with an initial tangential velocity. The particle decelerates due to a constant friction force once the particle has "settled" vertically. However, the particle unexpectedly experiences a negative tangential acceleration resulting in a negative tangential velocity, which does not appear physical. I would expect the effective friction force to be zero once the velocity becomes zero and hence there is no further acceleration of the particle (i.e. the particle will have "stuck" to the plane, so although there is a static friction force the velocity vector is 0).

test case:
test-friction.txt

velocity plot:

We would not expect to need to apply any tangential damping to the model to negate this behaviour. This will become complicated with the way we currently model "springs" and "dampers" since these are defined using separate SPRING definitions (e.g. see the test-friction.txt).

What do you think? Is it possible to change the SPRING definition such that stick/slip behaviour is modelled without defining tangential spring/damper values (ks/kn)?

Best regards,

Ben

[tkoziara]

Let me have a look at this example and get back; A simpler friction model is of course possible: one that applies friction always opposing slip velocity; In this case there is no need for tangential springs and dampers - but also there is now way to model frictional sticking and stick-slip transition; perhaps we can have a way of including both approaches and switching between them:) I'll have a look.

[tkoziara]

Ben - I have modified the model now so that kskn = 0 is allowed in the SPRING command (see updated docs). Practically the same effect can be obtained either by using kskn = 0 or by using kskn = very large number. See figure below.

There is an oscillation though when the particle comes to rest. When I print out friction force components it keeps switching between positive and negative every time step. You can see the pattern in the velocity here:

Why did you use two springs in the example? (the contact is undamped) When enabling damping in the contact spring and disabling the other spring I am getting awkward results:) So perhaps it is a good example for testing this. Thank you. Please, feel free to investigate from your end. I will also continue looking and let know when there is more. Do use the latest commit to test the above.

[bencannell]

Hi @tkoziara,

Thank you for looking into this and providing the figures. The behavior now looks more similar to what we expect.

I included two vertical springs in the model since this is how we model "contact" interactions in our large array models. One of the SPRING commands defines just a spring curve (the damping curve is set to [0.0, 0.0,...]), we then use a second SPRING command to represent a damper. I think this gives us more control over how the damping is applied and it is a direct translation of how the spring/damper pairs are modelled in the LS-DYNA models.

I'll let you know how I get on with the updated source code.

Ben

[bencannell]

Hi @tkoziara,

Just an initial observation, but the magnitude of the high frequency oscillations between +ve and -ve velocity appear to be very timestep dependent:

timestep = 0.5 * CRITICAL()

timestep = 0.1 * CRITICAL()

[tkoziara]

Yes, indeed, I also noticed that; Basically the friction force keeps switching from -f_tiny to +f_tiny while the magnitude of it is time-step dependent; So this seems to be a numerical artefact of some sort; I'll try to understand and eliminate it; I will look into creating a one spring model which reproduces that - if you were able to construct another one and paste it in - that might be helpful:)

[tkoziara]

Oh, look, I fixed a little bug (e0cdc1f), and with this example:
single-contact.py.gz
we can get these figures:



So the issue is now isolated with a simplest possible test:) I am wondering whether it will be possible to "fix" it or whether it is a modelling error pertinent to this type of modelling. I will share more when some conclusions come; (by the way, below I accidentally closed and reopened the issue)