Please Help

Trivolution

CP2Notes-ContactonCurvedSurface.pdf

Home >Engineering homework help >Civil Engineering homework help >Please Help

CEE 212--Dynamics

CP 2 Notes Contact on a curved surface

The Mechanics Project Arizona State University

Keith D. Hjelmstad

Contents. 1. The problem. 2. What is a plane? 3. Equation of motion. 4. Model of impact. 5. Practical modeling issues. 6. Finding the time of impact. 7. Updating for impacts in MATLAB 8. Multiple constraints 9. Summary.

CEE 212—Dynamics

-10 -5 0 5 10

-20

-15

-10

-5

CEE 212--Dynamics

CP 2 Notes The problem

The problem. The idea behind this project is to introduce contact problem by confining the motion of a projectile to be on one side of a surface. In essence, the particle will bounce off the walls of the surface. Between impacts, the motion of the particle is ordinary projectile motion, affected only by gravity and drag forces (but no geometrical constraints). To solve the contact problem we need to describe the surface (i.e., the location of the walls), establish a method to detect contact, and do an appropriate computation to resolve the impact problem. The geometries we will consider are pretty simple. We have surface defined as g(x,y) = 0, as shown in the sketch. The sketch shows an elliptical shaped surface. The particle bounces off the surface with a reduction in velocity using a model with a coefficient of restitution, which is a convenient way to deal with the energy lost in impact. What we will do in these notes is sort out the basic ideas of the contact problem by looking at what happens in a single contact with a single plane. The CP problem is then extend the idea curved surfaces.

A surface is defined by a

function

-10 -5 0 5 10

-20

-15

-10

-5

CEE 212--Dynamics

CP 2 Notes What is a plane?

What is a plane? A plane is a geometric object that embodies the idea of flatness. There are a few features of a plane that define it. First, and foremost, there is a unit vector n that is normal to the plane. What that means is that any vector that lies completely in the plane is perpendicular to the vector n. It should be evident that every plane parallel to our plane also has n as a normal vector. We can specify our plane by giving the location of one point in the plane, e.g., point A, which is located with position vector xA relative to the origin O. One way to specify a plane is to write the equation of the plane. In two dimensions (i.e., z = constant)

( , ) 0g x y y a x b   

where a and b are constants. If the constants are specified then the equation identifies a unique plane. The gradient of g is normal to the plane, i.e.,

g g

 

 n

The gradient of g can be computed (for the 2D case) as

1 2 1 2 g g

g a x y  

       

e e e e

One very important feature of the function

( , )g x y y a x b  

Is that it is positive for points on one side of the plane and negative for points on the other side. We will use this fact to detect contact.

O x

y ( , ) 0g x y 

( , ) 0g x y 

( , ) 0g x y 

plane

CEE 212--Dynamics

Equation of motion. We will divide all contact problems into three phases: Phase I—prior to contact, Phase II— during contact (which is very brief!), and Phase III—after contact. During contact (Phase II) the free body diagram is as shown at left. The equation of motion (balance of linear momentum) is

2W R m  e n x

We can multiply both sides of the equation by the projection P  InnT (see box at left) and integrate from 0 to 0+ to get

0 0 0

20 0 0

(0 ) (0 )

W dt R dt m dt

  

  

 

  

   

  P e P n P x 0 P x x



 

[ ] [ ] ( )

 

    

Pn I nn n I n nn n n n n n n n 0

The first term vanishes because the weight force is constant. The second term vanishes because Pn=0. The resulting equation is a statement of conservation of momentum. Note that it only applies to the part of the velocity that is in the plane and not to the part that is normal to the plane. Conservation of momentum tells us that the component of velocity in the plane does not change. What we do not know is how the velocity will change in the direction normal to the plane. Because we do not know the exact nature of the impact force R (i.e., how it varies with time during the very brief impact duration), we will need to make an assumption about the velocity in the normal direction.

CP 2 Notes Equation of motion (Phase II)

The matrix P I nnT is called a projection because it annihilates the n component of any vector. Observe its action on n itself:

So pre-multiplication by P clips off the n component of any vector and leaves alone the part of that vector in the plane normal to n.

R n 2W e

CEE 212--Dynamics

v  v n Pv

Model of impact. We know the velocity before impact, i.e., 0 , completely. What we want to find is the velocity after impact, i.e., 0 . We know we must leave the in-plane part the same and we need to find the normal part. So let us write the new velocity as

Where v is, as yet, unknown. Note that, by construction, we have guaranteed Pv+ = Pv, so no matter what value we get for v, we will still conserve momentum in the plane. Now we need a model. Let us assume that the normal component of the velocity after is in the opposite direction and is a fraction of the normal component of the velocity before impact, i.e.,

The quantity e is sometimes called the coefficient of restitution. The final velocity after impact should have this normal velocity and the previous tangential velocity. The coefficient of restitution would ideally be a sort of material constant that one could look up in a handbook. However, the local interactions during contact are very complex. Suffice it to say that this characterization is a good qualitative approach with the main value being to explore the effect of e between its two bounds.

CP 2 Notes Model of impact

  0 1v e e    n v

We can compute the energy lost during contact. Note that

      

  

1 2

2 1 2

1 2

221 2

E m

m e

  

  

  

 

   

 

   

v v

n v Pv Pv

v v

n v Pv Pv

 2212 (1 )lossE m e   n v Hence, the energy lost is

Note that energy lost depends upon the velocity of impact.

v v

CEE 212--Dynamics

CP 2 Notes Practical modeling issue

v v

Practical modeling issue. For the MATLAB code Contact_Plane.m (included as a template code for CP 2), the impact is modeled by simply testing the sign of the plane function g to see when it goes from positive (inside) to negative (outside). But, we need to look for is the exact condition that has

( ) 0 and ( ) 0old newg g x x

It should be theoretically possible to simply test the new position to see if g(xnew) < 0, because if you tested the previous step, and it was not negative, then that is the same as knowing that g(xold) > 0. But a practical problem lurks behind that logic. Imagine a case where the new position penetrates the plane and the coefficient of restitution is small enough that the new velocity is not adequate to get back to the positive side. You change the velocity and get ready for the next step. Now at the next step you still have g(xnew) < 0, but in this case g(xold) < 0, too. The algorithm would change the velocity and you would lose more energy. This is “getting stuck behind the glass.” It is not a physical contact issue, it is a numerical analysis issue. It is worthwhile to experiment with this case to convince yourself that reducing the time step will not always solve the problem. You really do need to check both states.

1nx

2nx

Numerical phenomenon of getting stuck under the glass ceiling. The restitution velocity is not enough to move xn+2 out, so algorithm computes red line reversal rather than green line advance. If you test both g(xold) > 0 and g(xnew) < 0, then this does not happen.

CEE 212--Dynamics

CP 2 Notes Finding the time of impact

Finding the time of impact. In a time-stepping algorithm we are not likely to hit the contact point on one of our discrete times. We can detect contact by testing g on all surfaces at each step. At the time where g goes from positive (inside) to negative (outside) then we need to make a contact computation. All we need is a test for positivity of the g function. Once we detect contact we can interpolate the position, velocity, acceleration, and time to give us the point of contact C on the surface. The fraction of the distance from the old point A (inside) to the new point B (outside) that just gets to the surface can be computed as

( ( )) 0old new oldg   x x x

Then we can linearly interpolate the new values

        

new old new

t t t t 

 

    

  

x x x

v v v

a a a

The linear interpolation is not exact, but should give good results.

In a time step in which contact is detected, we can simply adjust the values by ‘interpolation’ as indicated at right. These values then become the “new” values for the current time step. Note that when the impacts are very frequent (as in the later portions of the motion) the time advancement slows down a lot because the time step gets clipped so often.

oldx

newx

( ) 0oldg x

( ) 0newg x

contactx C



O x

y ( , ) 0g x y 

CEE 212--Dynamics

CP 2 Notes Updating for impacts in MATLAB

Updating for impact in MATLAB. We can compute the time, position, velocity, and acceleration of impact by solving g as

1nx

2nx

contactx

In MATLAB that looks like:

xold xnew

n



x x

( ( )) ( ) 0old new oldg g    x x x 

This is pretty easy for the linear function, but we must use Newton’s method (or some other iterative nonlinear equation solver) for curved surfaces. Newton’s method is

1 ( ) ( )

i i i

g g

  

   

 x

 

where x is xnew – xold (which is constant while we are finding ). We can start the iteration with z=0 since the value of z must be between zero and one. Sometimes this starting point will fail and it may be necessary to find a strategy to get convergence. Alternatively, we could use bisection or regula falsi (method of false position) to solve the nonlinear equations. The MATLAB function plane.m takes in the current value of the position of the particle x (which in the iteration is the most up-to-date value) and the parameters needed to compute the function g. It returns the value of the function and the gradient of the function at that position.

dx = xnew - xold; frac = 0; err = 1; it = 0; while (err > tol) && (it < itmax)

it = it + 1; x = xold + frac*dx; [g,dg] = plane(x,param,type); frac = frac - g/dot(dg,dx); err = abs(g);

end

CEE 212--Dynamics

CP 2 Notes Multiple constraints

1nx

1 0g 

2 0g 

1nx

2 1( )ng x

1 1( )ng x

Multiple constraints. If we have a corner created by multiple constraints it may happen that the particle passes out of the box and violates more than one of the constrains (as pictured in the sketch at left). When this happens we need some way to figure out which of the sides it would have contacted had we not passed right through it numerically. We can sort this question out by examining all of the constraints and picking the most violated one. Let us assume that we have arrived at location xn at time tn. At this point we check all of the constraints and find that

( ) 0 for all constraintsi ng ix

So we take another step. Now we have arrived at xn+1 and we check all of the constraints to find that we have violated two of them. To wit,

1 1 2 1( ) 0 ( ) 0n ng and g  x x

Both walls are candidates for the impact. We select the most violated one, which is wall 1, because

1 1 2 1( ) ( )n ng g x x

If more than one constraint is violated then we take the most violated constraint as the one to work with (and ignore the other one). All interpolated values are based on that one constraint.

CEE 212--Dynamics

CP 2 Notes summary

Summary. The phenomena of contact and impact are very important to a wide variety of engineering problems. In an earthquake, for example, two buildings that have been built close together will pound against each other and can cause lots of damage. Hence, pounding is an important design issue for tall buildings. Almost all dynamics loads begin with some sort of contact (note that in Statics we assume that the loads are just there and always have been there—they have to get there somehow…). This project explores the outcome of a single particle with multiple planes of constraint (the box). Keys to contact are (1) The detection algorithm (using what we call the gap

function, f). All we need to do is compute f at each step and look for when it goes from positive to negative.

(2) The interpolation of the state variables and time to exactly identify the state and time of contact. The linear interpolation of all three of the state variables is an approximation. (We could iterate on the time of contact and us the generalized trapezoidal rule to integrate from a to v to x.)

(3) The contact resolution computation using the model that respects conservation of momentum in the tangential direction but allows loss of energy in the normal direction.

If you think contact and impact are not important, just ask a dinosaur. Oh, wait, you can’t because impact wiped them out. Well, trust me, it’s important.