yade-dev team mailing list archive

[kan <nicessgg@xxxxxxxxx>]

Hi, Dear Bruno,

Thank you very much.

I have a typo in the last email.
Yes, the code is an increment pattern, the '+' is not input in here, but in
the code there is:
 Fs + = - Ks*double_ds;
vector_Fs=Fs* vector_Shear_Direction;
with debug model, I output the variables, so I have checked the  normal
directional vector, shear directional vector, and the force calculation for
Fn and Fs for each contact, several points are for sure:
1. the normal and shear directional vector are right;
2. the Fn and Fs are also right,

probably problem points may from
1. the sign error, but the possibility is low because I checked with debug
model and out put all the calculation steps for the first several steps, the
result looks resonable from the rotational direction;
2. I found when rotational speed is low, the shear relative speed is major
part is from the translational relative velocity's decomposed part along the
shear direction (that is to say, dot product with
shear_direction_unit_vector). since my Kn and Ks is chosen the same, so when
their is no normal relative velocity, the shear
delta_displacement ds=relative_Vs*dt, and the Fs+= Ks *ds, the amount is
huge. see the below calculation:
Assume Kn=Ks=1e9 N/m, particle radius is 1e-3 m, one particle is static,
another one is 20m/s along the shear direction, and both of them have no
initial rotation , density is 2.6e3 kg/m^3, then for 2D case--Disk with
thickness t=1 m, the mass is m=Pi*r*r*t*density= 8.16e-3 kg,  I=Inertia =
0.5*m*r*r= 4.08e-9 =4e-9, use dt=alpha*sqrt(m/k), the dt choose dt=1e-6,
then for the first iteration:
ds=vs*dt=20*1e-6=2e-5 m;           //vs has no rotational velocity part now
dFs=ds*Ks=(2e-5) * ( 1e9)=2e4 N; // (along the shear direction, so here use
the '+' )
Fs = Fs+dFs=0+2e4=2e4 N         // (along the shear direction, so here use
the '+' )
Moment=Fs*r=(2e4)*(1e-3) =20 N.m
rotational_acceleration =Moment/Inertia= 20/(4e-9) =5 e9 rad/sec
so after one iteration, the rotational speed is:
rotational_speed= rotational_speed + rotational_acceleration * dt = 0 +
(5e9)*(1e-6)= 5e3 rad/sec !!!!!! WRONG!!!!
here the damping force is not considered, or if we consider damping is 80%
of the resultal moment, the rotational_speed is still 1e3 rad/sec which is
still a large number .
Then what is the wrong point? too large Ks? Too large relative initial
speed(we need to know the translational speed may be increased to a large
value after long time simulation )? too large dt? or what else?

I checked with the CohesiveFrictional class in the file
'CohesiveFrictionalContactLaw.h, .cpp', in function 'void
CohesiveFrictionalContactLaw::action(MetaBody* ncb)', the calculate logic is
nearly the same, but why wo I get a wrong number?
Need your help. Thanks very much in advance.

//I CCed the email to yade-dev and yade-users, so some one who meets the
same problem could get help too, is it ok?

Kan



On Mon, Mar 2, 2009 at 7:07 AM, Bruno Chareyre
<bruno.chareyre@xxxxxxxxxxx>wrote:

>  Hello
>
> Is the pseudo code below for developping a new contact law? It is probably
> useless since there is one law doing what you need already : the
> "CohesiveFrictional" class.
>
> There are probably two mistakes in your code :
> 1/ For sure : what you define as Fs= - Ks*double_ds is in fact a dFs (an
> increment), not the total force.
> If you check the current contact laws, you will see something like Fs = Fs
> - Ks*double_ds. With your equations, Fs is in fact a viscous force,
> proportional to the current relative velocity.
>
> 2/ Not sure : the fact that rotational velocity is constantly increasing
> may be due to a sign error somewhere (e.g. in the normal orientation), or
> perhaps an inapropriate timestep (too large). I don't feel like the first
> mistake could produce this result.
>
> Bruno
>
>
> Nice Man a écrit :
>
>> Hi,Dear all.
>>  I would like to use cohesive contact model to simulate a rock, the small
>> rock particle is modeled by cementing them and the force assumed by the
>> cementing is cohesive force--which can assume normal and shear force.  but I
>> meet problem in moment calculation and the rotational speed calculation:
>> just for 2D case ( 'x' is cross(), Normal_directior_unit is from P1 center
>> to P2 center ):
>> vector_relative_velocity=(Particle1.v+Particle1.omega x
>> (Particle1.r*Normal_directior_unit) ) -(Particle2.v+Particle2.omega x
>> (Particle2.r * -Normal_directior_unit) ) ;
>> vector_relative_shear_velocity=(vector_relative_velocity *
>> vector_Shear_Direction) * vector_Shear_Direction;
>> double_relative_shear_velocity=(vector_relative_velocity *
>> vector_Shear_Direction)
>> double_ds=double_relative_shear_velocity*dt;
>>  Fs= - Ks*double_ds;
>> vector_Fs=Fs* vector_Shear_Direction;
>>  Now moment is :
>> delta_Moment= (Particle1.r * Normal_directior_unit) x vector_Fs; //(Fn is
>> not calculated here,and also since Fn pass through the center of the mass,
>> so it has no effect on the rotation)
>>  Total_moment += delta_Moment;
>>  so
>> Rotational_velocity += ( Total_moment / Inertia)*dt;
>>  My problem met is , after several iterations,
>> my rotational velocity reachs 1e4 rad/sec, or even more!!! this is totally
>> incorrect.
>> So what is the calculation problem in this cohesive condition? ?  or can I
>> do the moment calculation and rotational speed calculation like this since
>> there is cementing cohension between particles and the cohension will strict
>> the free to rotate of the particle now? In the CohesiveFrictionalContactLaw,
>> I see the moment calculation.
>>  need your help. Thanks in advance.
>>  Kan
>> ------------------------------------------------------------------------
>>
>> _______________________________________________
>> Mailing list: https://launchpad.net/~yade-dev
>> Post to     : yade-dev@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~yade-dev
>> More help   : https://help.launchpad.net/ListHelp
>>
>>
>
>
> --
>
> _______________
> Chareyre Bruno
> Maitre de conference
>
> Grenoble INP
> Laboratoire 3SR - bureau E145
> BP 53 - 38041, Grenoble cedex 9 - France
> Tél : 33 4 56 52 86 21
> Fax : 33 4 76 82 70 43
> ________________
>
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-dev
> Post to     : yade-dev@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~yade-dev
> More help   : https://help.launchpad.net/ListHelp
>