W elcome to
Fintronic USA

redball.gif (326 bytes)About Fintronic USA

redball.gif (326 bytes)Main announcements

redball.gif (326 bytes)What's new at Fintronic

redball.gif (326 bytes)What our customers are saying...

redball.gif (326 bytes)Support for SystemC

redball.gif (326 bytes)Support for Verilog 2001

redball.gif (326 bytes)Third party tools integrated with FinSim(Specman, Denali, Debussy, Undertow, Vera, VirSim, HDL Score, Comet, Meteor, NelSim, Pivot, DeskPOD, @Designer 2.1)

home.htmlProductsSupportEvaluateContact

Differential Eq. with variable coefficient in Verilog: model of the force of a tennis ball hitting a wall


This example works on FinSim 10_01_31 and subsequent versions.

This example models the force of a tennis ball hitting a wall, based on ideas presented by S J Haake, M J Carre, R Kirk, and S R Goodwill at the University of Sheffield, Dept. of Mechanical Engineering.. The main model is that of a spring-mass system governed by the formula: m*y(2) + c*y(1) + k*y = 0, with y = 0 and y(1) = speed of ball hitting the wall. There are three forces that push against the wall:

1) the spring force, spring_force is due to the ball being compressed, with spring_force = k*y.

2) the damper force, damper_force is due to the viscuosity of the ball with damper_force = c*y(1)

3) the flux force, flux_force is due to the portion of the ball that loses it's speed during a small time slice. This flux_force is m0 * speed / delta_time, where m0 = m*s1/s2, where s1 is approx. the surface of a cillinder having as base the flattened portion of the ball and as hight the current speed y0 * delta_time and s2 is the surface of the ball, which makes

flux_force = y0*y0*m* $VpSqrt(2*r*x0-x0*x0)/(2*r*r), where y0 is the speed at a particular time, m is the total mass of the ball, r is the radious of the ball, and x0 is the position of the center of the ball with respect to its original position when the impact began.

The algorithm solves the problem for speeds of 10m/s and 30m/s in a loop. In an inner loop the algorithm adjusts the non-linear coefficients and solves the equation for a number of steps with constant coeficients. It then performs again the same steps for the new value of the loop variable (named slice), which automatically adjusts which portion of the data is provided to the differential eq. solver via the "view..as" mechanism. Note that by using the "view..as" mechanism no data transfer is needed during the computation, the data being computed "in place".

The results are displayed via the $PrintM task as well as via the $VpPtPlot task. The latter produces the results shown here 

module top;
parameter real alpha = 1.65;
parameter real r_total_time = 0.005;/* seconds */
parameter nr_pts_per_ct_coef = 10;
parameter nr_slices_per_sec = 10000;
parameter nr_pts_per_sec = nr_pts_per_ct_coef * nr_slices_per_sec;
parameter integer Size = r_total_time*nr_pts_per_sec;
parameter real r_size = Size;
parameter real h = 2*r_total_time/r_size;/* set double of sampling period */
parameter real m = 0.057;/* kg */
parameter real A = 16000000;/* N/m**2 */
parameter real k0= 21000; /*N/m*/
parameter real B = 3500; /* Ns/m */
parameter real r = 0.032; /*m (radius of ball)*/
parameter order = 2;
parameter nrEq = 1;
localparam integer nr_slices = nr_slices_per_sec * r_total_time;

real Fe[0:0][0 : Size], x[0:0][0 : Size], y[0:0][0 : Size];
reg [0 : 3199] ressymb[0 : 0];
real coef[0 : 0][0 : 2];

view real x_ct [0:0][0:nr_pts_per_ct_coef] 
    as x[0][$I2+slice*nr_pts_per_ct_coef];
view real Fe_ct [0:0][0:nr_pts_per_ct_coef] 
    as Fe[0][$I2+slice*nr_pts_per_ct_coef];
view real y_ct [0:0][0:nr_pts_per_ct_coef] 
    as y[0][$I2+slice*nr_pts_per_ct_coef];

real ar_f_spring [0:nr_slices];
real ar_f_damper [0:nr_slices];
real ar_f_flux [0:nr_slices];
real ar_f_total [0:2][0:nr_slices];

integer i, j, isZero, total_time;
integer slice = 0;
real k, c, x0, y0, v0;

initial begin
  #1;
  coef[0][0] = m;
  $InitM(Fe, 0.0);
  for (j = 0; j < 2; j = j + 1)
  begin
    /* find total force for two speeds: v0 = 10m/s and v0 = 30m/s */
    v0 = (j == 0) ? 10 : 30;

    x[0][0] = 0;
    y[0][0] = v0;
    isZero = 0;/* kludge for nulifying results for when the model does not apply*/
    for (slice = 0; slice <= nr_slices; slice = slice + 1)
    begin
      x0 = (x_ct[0][0] > 0) ? x_ct[0][0] : -x_ct[0][0];
      y0 = y_ct[0][0];
      c = 4.0*B*x0*(2.0*r-x0);
      k = (slice < 2) ? 120000 : (k0 + A*(x0**alpha));

      coef[0][1] = c;
      coef[0][2] = k;
      ar_f_spring[slice]= x0*k;
      ar_f_damper[slice] = y0*c;
      if ((y0 > 0) && (x0 > 0)) begin
           ar_f_flux[slice] = y0*y0*m*
                              $VpSqrt(2*r*x0-x0*x0)/(2*r*r);
      end
      else  ar_f_flux[slice] = 0;

      ar_f_total[j][slice] = ar_f_spring[slice] + 
                             ar_f_damper[slice] +
                             ar_f_flux[slice];
      ar_f_total[j][slice] = (isZero) ? 0: ((ar_f_total[j][slice]>0)?ar_f_total[j][slice]:0);
      if ((slice != 0) && (ar_f_total[j][slice] < 15 
           /* after this value the spring model no longer applies*/ )) begin
        isZero = 1;
      end
      /* call dif eq solver. The array x will contain the solution and the
       array y will contain the first derivative of x.
       Note that the initial conditions have been already placed in x[0]
       and y[0]. 
       */
      $VpLODE(order, nrEq, h, nr_pts_per_ct_coef+1, 
              x_ct, coef, Fe_ct, y_ct, ressymb);
    end
    ar_f_total[j][nr_slices] = 0.0;
  end
  $PrintM(ar_f_total,"%e");
  total_time = r_total_time*1000;
  $Flot("ball_force.html", 2, h/2, 
           "Tennisball_(0.057kg,0.032m)_Force_pushing_the_wall)", 
           "Time (ms)", "Total Force(N)", 0, nr_slices, ar_f_total,
           "10m/s", "30m/s");
 end
endmodule

© Copyright 1999-2021, Fintronic USA, Inc.   All rights reserved.