Calculating change in velocity with resistance

[kennethgilpin]

Hi everyone,

I'm a doctor who is fair at maths - haven't done physics for a while. I'm trying to model some heart physiology which I can describe as a mechanics problem:

I have a particle, mass m, with a constant force acting on it, f.

It also has resistance R, where the resistance is equal to the velocity of the particle squared times a constant k:
R = kv²

my equation of motion is:

m.dv/dt = F = f - kv²

I want to find the velocity after a given time. The are formula on the internet which describe a solution to this problem (mostly in the context of a falling body accelerating to wind resistance):

http://keisan.casio.com/has10/SpecE...sistance (velocity and distance)/default.xml"

However they all assume that starting velocity is zero. I need a solution for any starting velocity. so:

t₀ = time 0.
v₀ = velocity at time 0.
f = constant force on particle
m = mass
R = resistance = kv² as above, where k = constant, v = velocity

what would be the velocity at time t?

I know this is difficult. I'm not sure if this is really a maths problem.

kenneth

 

[Filip Larsen]

In short, you have two cases: either the speed will start below the terminal speed and increase or it will start above terminal speed and slow down. The terminal speed is where the acceleration is zero,

(1) v_t = \sqrt{f/k}

so for v_0 < v_t the solution for v is

(2) v = v_t \tanh(\frac{k v_t}{m} t + \textrm{arctanh}( \frac{v_0}{v_t} ) )

where \tanh and \textrm{arctanh} are hyperbolic tangent and its inverse, and for v_0 > v_t it is

(3) v = v_t \coth(\frac{k v_t}{m} t + \textrm{arccoth}( \frac{v_0}{v_t} ) )

where \coth and \textrm{arccoth} are hyperbolic cotangent and its inverse.

 

[Filip Larsen]

By the way, the distance traveled will in both cases end up being expressed by

(4) s= s_0 + \frac{m}{k}\log\left[\frac{v_0}{v_t}\sinh(\lambda t)+\cosh(\lambda t)\right]

where I've introduced the constant

(5) \lambda = \frac{k v_t}{m} = \frac{\sqrt{fk}}{m}.

 

[kennethgilpin]

thank you for your quick reply. The function seems to work well in with a few test examples I've just put in.

the problem is for my use is with deceleration due to a change in direction of the force, f.

if I have a positive velocity, v₀, but a negative force, f, then function (2) does not provide a v less than then velocity v₀.

clearly if the terminal speed is related to the root of the force then it cannot be negative. If I use the absolute force and multiple the terminal speed by the direction of the force, i.e.

v_t = |f|/f * \sqrt{}|f|/k

then the terminal speed has the same sign as the force direction, which makes sense. However, when this number is put into function (2) the subsequent value, v is still greater than v₀.

any ideas?

k.

 

[Filip Larsen]

Ok, I assumed f was positive. If f is negative, you would normally write the differential equation as

(6) m\dot{v} = -f - k v^2

with f positive. This equation has the solution

(7) v = v_t \tan(\arctan(v_0/v_t)-\lambda t)

for initial speed v_0 and the time in the range 0 \le t \le t_1, where

(8) t_1 = \arctan(v_0/v_t)/\lambda

is the time until the particle has reached zero speed. After that time the movement of the particle is governed by (2) with initial speed zero. The position for (7) in the same time range is

(9) s = s_0 + \frac{m}{k}\log\left[\frac{v_0}{v_t}\sin(\lambda t)+cos(\lambda t)\right]

 

[kennethgilpin]

thank you all of the formula work wonderfully.

You haven't defined v in the case of f being negative, and also v₀ > v_t. I assumed it might be:

<br /> v = v_t \cot(\cot^{-1}(v_0/v_t)-\lambda t) <br />

(can't find the Latex for arccot!)

but this doesn't work out, and when I graph the solution the gradient is positive (i.e. the solution, v, is larger than v₀). Could you help me out for the last time?

increasingly gratefully,

kenneth

 

[Filip Larsen]

You haven't defined v in the case of f being negative, and also v₀ > v_t.

I didn't explain it properly. For a negative force (i.e. here a force that oppose the initial direction of speed), you can use (7) no matter what the initial speed is. However, (7) can only be used until the speed reaches zero (which it does after the time shown in (8)) since for negative speeds the fixed minus in front of the drag force will make drag force act in same direction as the speed, which is not physical.

So, when the speed reach zero using (7), you will have to switch to using (2) (which handles the case where speed and force point in same direction for low initial speeds) with initial speed zero and you will now have to consider the speed calculated from (2) as being in the opposite direction as it was when using (7).

For example, let's say your particle moves so that positive speed is to the right and with the force pointing left. You use (7) until v = 0 and then use (2) with v₀ = 0 and where you now consider positive speeds pointing left (and the force still pointing left). If you are interested in the distance traveled too you will also have to carry over the distance from (9) as initial distance in (4) when switching equations and, like with the speed, consider the distance to now increase in the other direction.

I hope its a bit more clear now otherwise just ask.

 

[kennethgilpin]

Filip thank you a thousand times for your help with this problem.

I just wanted to tie this thread off in case someone else comes here looking for a solution. You can see that I've used some of the formula above but haven't written in the need to decelerate (or accelerate) the particle to zero when changing direction. I think this is because I'm using these formulae in a loop in a computer program which considers small intervals of time. Filip Larson might want to comment on the code. This snippet is in an ECMA script language:

var dt:Number = 0.005; //seconds per interval, so 200 intervals = 1sec
var i:int=0;
var f:Number = 10;  //force
var v:Number = 0;   //velocity - starting from stationary here
var k:Number = 0.1; //constant
var m:Number = 5;   //mass

//determine terminal velocity
var Vterm:Number = Math.pow(Math.abs(f/k),0.5);

//introduced constant
var gam:Number = Math.pow(Math.abs(f * k),0.5) / m;

for (i; i < 600; i++)	{
//i.e. 600 / 0.005 = 3 secs

//is velocity above terminal speed?
if (Math.abs(v) < Math.abs(Vterm))	{
   if (dP > 0)	{	//accelerating
	v= Vterm * advMath.tanh(gam * dt + advMath.atanh(v/Vterm));
   }
   else	{	//decelerating
	v= Vterm * Math.tan(Math.atan(v/Vterm) - gam * dt);
   }						
}
else	{  //is velocity below terminal speed
   if (dP > 0)	{	//accelerating
	v= Vterm * advMath.coth(gam * dt + advMath.acoth(v/Vterm));
   }
   else	{	//decelerating
	var newVel:Number = Vterm * advMath.cot(advMath.acot(v/Vterm) - gam * dt);
	v +=  v - newVel;
        //this is clearly a fudge
   }
}	



}

where advMath.coth etc. are custom functions for the hyperbolic cotangents etc. See wikipedia for the formulae.

hope this helps someone. I'm a big fan of this forum

kenneth

 

[Filip Larsen]

I just wanted to tie this thread off in case someone else comes here looking for a solution. You can see that I've used some of the formula above but haven't written in the need to decelerate (or accelerate) the particle to zero when changing direction. I think this is because I'm using these formulae in a loop in a computer program which considers small intervals of time. Filip Larson might want to comment on the code.

This snippet is in an ECMA script language:

var dt:Number = 0.005; //seconds per interval, so 200 intervals = 1sec
var i:int=0;
var f:Number = 10;  //force
var v:Number = 0;   //velocity - starting from stationary here
var k:Number = 0.1; //constant
var m:Number = 5;   //mass

//determine terminal velocity
var Vterm:Number = Math.pow(Math.abs(f/k),0.5);

//introduced constant
var gam:Number = Math.pow(Math.abs(f * k),0.5) / m;

for (i; i < 600; i++)	{
//i.e. 600 / 0.005 = 3 secs

//is velocity above terminal speed?
if (Math.abs(v) < Math.abs(Vterm))	{
   if (dP > 0)	{	//accelerating
	v= Vterm * advMath.tanh(gam * dt + advMath.atanh(v/Vterm));
   }
   else	{	//decelerating
	v= Vterm * Math.tan(Math.atan(v/Vterm) - gam * dt);
   }						
}
else	{  //is velocity below terminal speed
   if (dP > 0)	{	//accelerating
	v= Vterm * advMath.coth(gam * dt + advMath.acoth(v/Vterm));
   }
   else	{	//decelerating
	var newVel:Number = Vterm * advMath.cot(advMath.acot(v/Vterm) - gam * dt);
	v +=  v - newVel;
        //this is clearly a fudge
   }
}	
}

where advMath.coth etc. are custom functions for the hyperbolic cotangents etc. See wikipedia for the formulae.

Some quick comments:

• If you want to tabulate the speed of a particle, you should probably calculate each speed value from the original solution using total time. In the code above, the analytical solutions are used to advance the speed from the speed at time t to the speed at time t+dt and then this speed is used as the initial speed for the next time step. Doing that will make the solution "drift" due to truncation and rounding error in the calculations. A numerically more stable technique is to calulate t = t₀+i*dt and v = v_solution(t), using the appropriate solution depending on conditions.
  
• The particle can be in four different regimes based on its initial speed and direction relative to the force. Classifying the regimes based on what type of function used for the speed solution we have that regime 1 corresponds to the tan solution, regime 2 to the tanh solution, regime 3 to the terminal speed solution, and regime 4 to the coth solution. A particle in regime 2, 3 and 4 will always stay the same regime, while a particle in regime 1 will eventually migrate to regime 2.
  
• It is probably easier to keep the force always positive and use the sign of the speed to indicate if it is opposite the force or not, for instance using negative speed when the particle is moving opposite the force.

If I modify your code to incorporate the above points, I get the following UNTESTED code:

// parameters
var dt:Number = 0.005; //seconds per interval, so 200 intervals = 1sec
var f:Number = 10;  //force - always positive
var v0:Number = 0;  //initial speed  - positive when same direction as f
var k:Number = 0.1; //constant - always positive
var m:Number = 5;   //mass - always positive
var t0:Number = 0;  // initial time
// state variables
var v:Number = v0;  //speed - positive when same direction as f
var t:Number = 0;   // time
var i:int=0;        // step number

//determine terminal velocity
var vt:Number = Math.pow(f/k, 0.5);

//introduced constant
var gam:Number = Math.pow(f*k, 0.5) / m;

//time parameter used in the four regimes
var tz:Number = t0;
if (v0 < 0) {   // regime 1 - speed opposing force
  tz = Math.atan(-v0/vt)/gam;
} else if (v0 < vt) { // regime 2 - speed non-negative and below vt
  tz = advMath.atanh(v0/vt)/gam;
} else if (v0 == vt) { // regime 3 - speed equal vt
  // nothing special
} else { // regime 4 - speed above vt
  tz = advMath.acoth(v0/vt)/gam;
}

for (i; i < 600; i++) { //i.e. 600 / 0.005 = 3 secs
  t = t0 + i*dt;
  if (v0 < 0) { // regime 1
    if (t < tz) { // still in regime 1?
      v = -vt * Math.tan(gam * (tz-t));
    } else { // swith to regime 2 solution
      v0 = 0; t0 = tz; tz = 0; 
      t = t0 + i*dt;
    }
  }
  if (v0 >= 0) {
    if (v0 < vt) { // regime 2
      v = vt * advMath.tanh(gam * (t+tz));
    } else if (v0 == vt) { // regime 3
      v = vt;
    } else { // regime 4
      v = vt * advMath.coth(gam * (t+tz));
    }
  }
  // print or use t and v here 
}

 

[kennethgilpin]

thank you for your reply. I've been using a mixture of your code and mine and it's been working well. I think that there is another regime and i wanted your opinion/help for it:

where, speed is above vt _and_ speed opposes the force.

once again you would still need a 'tz' time parameter for when/if the speed changes direction.

kenneth

 

[Filip Larsen]

There should be no need for more regimes. The tan solution (regime 1) is valid for all speeds in opposite direction of the force. The terminal velocity does not have direct physical meaning as a terminal velocity in this regime since the terminal is relevant only for regime 2, 3 and 4. The time parameter tz in regime 1 signifies the time it takes for the particle to change speed from -v0 to 0 and is theoretically valid for arbitrarily negative speeds.

 

[Filip Larsen]

I now notice that the code I included in #9 is wrong regarding switching from regime 1 to 2. It should be

    ...
    } else { // switch to regime 2 solution
      v0 = 0; tz = -tz; 
    }
    ...

With this change t0 is always zero and could be left out.