Message Boards Message Boards

Anonymous User

NDSolve Problem for Solving Nonlinear Sytem of 3 Ball Chain "Hertz Model"

Anonymous User
Posted 6 years ago
4 Replies
0 Total Likes

Hello There, I've trying to solve 3 nonlinear equation of 3 ball chain connected with springs using Hertz model. Simply 3 balls connect with two springs. I was successful in obtaining the differential equations with no problems. I am using steel balls so you might find parameters for material properties such as Elastic modulus, Poisson ratio, radii of the balls...etc. The main problem is that when using the NDSolve function I get that the system is either singular or stiff. I tried many different methods such StiffnessSwithing, and ExplicitEuler...etc but non of them work especially when I try to solve for larger amount of time such t = 10 seconds. When solving for very small time scale such as micro second that works but definitely not correct. My initial conditions are subject to change with no problem. I would really appreciate your help on this. I am pasting the code in here and I am also going to attach it for you in this post. I did comment on my code so you can understand it with no problems.

(*values of masses,Poisson ratios,modulus of elasticity, Radii of the balls*) 
    values = {m1 -> 1, m2 -> 2, m3 -> 3, E1 -> 210 10^9, E2 -> 210 10^9, 
       E3 -> 210 10^9, \[Nu]1 -> 0.3, \[Nu]2 -> 0.3, \[Nu]3 -> 0.3, 
       R1 -> 6.98*0.01/2, R2 -> 6.98*0.01/2, R3 -> 6.98*0.01/2};
    (*Define the displacement vector for each ball*)
    x[t_] := {x1[t], x2[t], x3[t]};
    (*Assumed initial conditions*)
    IC = {x1[0] == 0.0000001, x2[0] == 0, x3[0] == 0, x1'[0] == 1, 
       x2'[0] == 0, x3'[0] == 0};
    k1 = (1 - \[Nu]1^2)/(\[Pi] E1) /. values;
    k2 = (1 - \[Nu]2^2)/(\[Pi] E2) /. values;
    k3 = (1 - \[Nu]3^2)/(\[Pi] E3) /. values;
    (* Contact Stiffness *)
    K12 = 4/(3 \[Pi] (k1 + k2)) Sqrt[(R1 R2)/(R1 + R2)] /. values;
    K23 = 4/(3 \[Pi] (k2 + k3)) Sqrt[(R2 R3)/(R2 + R3)] /. values;
    (*Elastic deformation at contact points between the balls *)
    \[Delta]12 = (x1[t] - x2[t]) + (R1 + R2) /. values;
    \[Delta]23 = (x2[t] - x3[t]) + (R2 + R3) /. values;
    (* Hertz Equations: Contact forces between the balls *)
    F12 = K12 \[Delta]12^1.5;
    F23 = K23 \[Delta]23^1.5;
    (*Differential equations for each ball*)
    EQ1 = m1 x1''[t] == -F12;
    EQ2 = m2 x2''[t] == F12 - F23;
    EQ3 = m3 x3''[t] == F23;
    (*Put them all together *)
    EOM = {m1 x1''[t] == -F12, m2 x2''[t] == F12 - F23, 
       m3 x3''[t] == F23} /. values
    (*Solve the differential equations *)
    eqn = NDSolve[{EOM, IC}, {x1[t], x2[t], x3[t], x1'[t], x2'[t], 
         x3'[t]}, {t, 0, 10}, Method -> "StiffnessSwitching"] // Flatten //
    (* Plot displacement and velocity responses of the system *)
    Plot[Evaluate[x1[t] /. eqn], {t, 0, 10}]
    Plot[Evaluate[x2[t] /. eqn], {t, 0, 10}]
    Plot[Evaluate[x3[t] /. eqn], {t, 0, 10}]
POSTED BY: Anonymous User
4 Replies
Posted 6 years ago

If I suspect there might be numerical accuracy problems and thus I turn all your decimal floating point values into exact fractions, include the options , WorkingPrecision->64, AccuracyGoal->16 inside your NDSolve, restrict your plots to {t, 0, 0.000839} because that is where NDSolve says things blow up, include PlotRange -> All in each plot to make it more likely that I see all the behavior and separately plot x1, x1', x2, x2', x3, x3' then I see your derivatives smoothly grow from zero to as much as 120000 while t grows from 0 to to 0.000839. Your derivatives are about 20 billion times the size of your positions.

If you look at your system should those derivatives be that huge that quickly? Or is there perhaps an error anywhere in your system? Even using hundreds of digits of working precision doesn't change the behavior.

Anonymous User
Anonymous User
Posted 6 years ago

Dear Bill. Yes I was able to make work after I restricted the time to very small number. I tried your suggestion and it work too. Thank you. The only problem causing the code not work is the nonlinear power of 3/2 in the equations F12 and F23. If you set 3/2 to 1 "cancel the nonlinearity" in the system everything works perfectly but again it's not going to be able to give the results I am looking for since I linearized my system. I am not sure if there is any other way of solving since I have been trying to for a few days. Thank you again.

POSTED BY: Anonymous User
Posted 6 years ago

What happens if you replace all your decimal approximations with exact fractions so there will be no round-off error and you write your own simple Euler numerical solver? Your equations give explicit expressions for your second derivatives. It looks like you can get linear approximations of your first derivatives from those and then get linear approximations of your functions from those. Iterate that over each small step and see how it behaves.

I understand many users want a magic "do what I want" button. I realize there must be a better way, but if you can't understand and successfully control the inner works of any of the black box methods then writing a very simple numerical solver that you could understand might get you an answer today and that might be enough to get you started. That only seems feasible because the form of your system and initial conditions are simple enough to perhaps allow this. I've written a dozen lines to try to do that but I'm still trying to convince myself it is correct.

Anonymous User
Anonymous User
Posted 6 years ago

Dear Bill, I did what you suggested and indeed it works for larger time scale except for the plots they go up to 0.5 seconds which is definitely better than before. I turned my values to exact fraction that works perfectly. I believe that work is so far correct if you are still not convinced google the following journal paper with the following title "Shock Absorption Using Linear Particle Chains With Multiple Impacts" and check section 2.3 of the Hertz model. Do you suggest that there is any other way of solving the problem for larger amount of time ? What would you suggest in terms if coding ? Like making a numerical solver that solves this particular problem! Thank you so much you've been very helpful!

POSTED BY: Anonymous User
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract