Anonymous User

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

Anonymous User
Posted 6 years ago
5635 Views
|
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 // Chop (* 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}]  Attachments:
4 Replies
Sort By:
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 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.