Force, Moment, and Torque

Note

You can download this example as a Python script: loads.py or Jupyter Notebook: loads.ipynb.

import sympy as sm
import sympy.physics.mechanics as me
me.init_vprinting(use_latex='mathjax')

class ReferenceFrame(me.ReferenceFrame):

    def __init__(self, *args, **kwargs):

        kwargs.pop('latexs', None)

        lab = args[0].lower()
        tex = r'\hat{{{}}}_{}'

        super(ReferenceFrame, self).__init__(*args,
                                             latexs=(tex.format(lab, 'x'),
                                                     tex.format(lab, 'y'),
                                                     tex.format(lab, 'z')),
                                             **kwargs)
me.ReferenceFrame = ReferenceFrame

Learning Objectives

After completing this chapter readers will be able to:

• Calculate the resultant force acting on a set of particles or bodies.

• Calculate the moment of a resultant about a point.

• Find the equivalent couple of a torque and resultant to simplify forces applied to bodies.

• Determine if a force is noncontributing or not.

• Define sign conventions for forces.

• Define forces for gravity, springs, friction, air drag, and collisions.

Force

A force is an abstraction we use to describe something that causes mass to move (e.g. accelerate from a stationary state). There are four fundamental forces of nature of which all other forces can be derived. We describe the resulting effect of these fundamental forces in this text. Moments and torques arise from forces and are useful in describing what causes distributed mass rotation. Forces, moments, and torques have magnitude and direction and thus we use vectors to describe them mathematically.

Bound and Free Vectors

Vectors may have a line of action. A line of action is parallel to the vector and passes through a particular point. If a vector has a line of action, it is said to be bound to its line of action. If a vector is not bound to a line of action it is said to be free.

Angular velocity is an example of a free vector. It has a direction, sense, and magnitude, but is not associated with any line of action. For example, if a disc can rotate around a fixed point, you can place a body anywhere on this disc and the body will always have the same angular velocity. A force vector, on the other hand, is bound. If a force is applied to a rigid body, we must know where on the body it is applied to resolve the force’s effect. A force vector acting on rigid body \(B\) at point \(P\) has a line of action through \(P\) and parallel to the force vector.

Fig. 34 Depicition of bound a) and free b) vectors.

Moment

If a vector is a bound vector, then we can define its moment about a point. The moment \(\bar{M}\) of bound vector \(\bar{v}\) about point \(P\) is a vector itself and is defined as ([Kane1985], pg. 90):

(158)\[\bar{M} := \bar{r}^{L/P} \times \bar{v}\]

\(\bar{r}^{L/P}\) is a position vector from \(P\) to any point \(L_i\) on the line of action \(L\) of \(\bar{v}\). The cross product definition ensures that the the moment does not depend on the choice of the point on the line.

Fig. 35 \(\bar{v}\) is bound to a line \(L\). The moment can be calculated based on a position vector from \(P\) to any point on the line, for example \(L_1,L_2\) or \(L_3\) as shown.

A moment can be the result of a set of vectors. The resultant of a set \(S\) of vectors \(\bar{v}_1,\ldots,\bar{v}_\nu\) is defined as:

(159)\[\bar{R}^{S} := \sum_{i=1}^{\nu} \bar{v}_i\]

If each vector in the resultant is bound, the sum of the moments due to each vector about \(P\) is called the moment of \(\bar{R}^{S}\) about \(P\). This summation can be written as:

(160)\[\bar{M}^{S/P} = \sum_{i=1}^{\nu} \bar{r}^{L_i/P} \times \bar{v}_i\]

where \(L_i\) is a point on the line of action of the associated \(\bar{v}_i\).

The moment of the set of bound vectors \(S\) about one point \(P\) is related to the moment of the set about another point \(Q\) by ([Kane1985], pg. 91):

(161)\[\bar{M}^{S/P} = \bar{M}^{S/Q} + \bar{r}^{P/Q} \times \bar{R}^{S/Q}\]

where \(\bar{R}^{S/Q}\) is the resultant of the set \(S\) bound to a line of action through point \(Q\).

For example, take the set \(S\) of two bound vectors \(\bar{F}_1\) and \(\bar{F}_2\) bound to lines of action through points \(P_1\) and \(P_2\), respectively. Below I’ve given the vectors some arbitrary direction and magnitude.

N = me.ReferenceFrame('N')

F1 = 2*N.x + 3*N.y
F2 = -4*N.x + 5*N.y

r_O_P1 = 2*N.x
r_O_P2 = 3*N.x

\(\bar{M}^{S/P}\) can be calculated directly using Eq. (160):

r_O_P = -5*N.x

M_S_P = me.cross(r_O_P1 - r_O_P, F1) + me.cross(r_O_P2 - r_O_P, F2)
M_S_P

\[\displaystyle 61\hat{n}_z\]

Or if \(\bar{M}^{S/Q}\) is known, as well as \(\bar{r}^{P/Q}\), then the Eq. (161) could be used:

r_O_Q = 5*N.y
M_S_Q = me.cross(r_O_P1 - r_O_Q, F1) + me.cross(r_O_P2 - r_O_Q, F2)

M_S_P = M_S_Q + me.cross(r_O_Q - r_O_P, F1 + F2)
M_S_P

\[\displaystyle 61\hat{n}_z\]

Couple

A set \(S\) of bound vectors with a resultant equal to zero is called a couple. A couple can have as many vectors as desired or needed with a minimum number being two, such that \(\bar{R}^{S}=0\). A couple composed of two vectors is called a simple couple. Fig. 36 shows a few examples of couples.

Fig. 36 Three couples: a) simple couple, b) & c) couples made up of multiple forces

The torque of a couple, \(\bar{T}\), is the moment of the couple about a point. Because the resultant of a couple is zero it follows from (161), the torque of a couple is the same about all points. The torque, being a moment, is also a vector.

Equivalence & Replacement

Two sets of bound vectors are equivalent when they have these two properties:

1. equal resultants

2. equal moments about any point

If 1. and 2. are true, the sets are said to be replacements of each other. Couples that have equal torques are equivalent, because the resultants are zero and moments about any point are equal to the torque.

Given a set of bound vectors \(S\) and a set of bound vectors that consist of a torque of a couple \(\bar{T}\) and vector \(\bar{v}\) bound to an arbitrary point \(P\) it is a necessary and sufficient condition that the second set is a replacement of the first if ([Kane1985], pg. 95):

(162)\[\begin{split}\bar{T} = \bar{M}^{S/P} \\ \bar{v} = \bar{R}^{S/P}\end{split}\]

This means that every set of bound vectors can be replaced by an equivalent torque of a couple and a single bound vector that is the resultant of the replaced set. This replacement simplifies the description of forces acting on bodies.

Take for example the birds eye view of a four wheeled car which has front steering and motors at each wheel allowing for precise control of the propulsion forces at each wheel. A diagram of the forces acting at each wheel is shown in Fig. 37.

Fig. 37 Set \(S\) of forces acting at each tire can be replaced with a resultant and a torque at a specified point, in this case \(B_o\).

In SymPy Mechanics, first define the symbols:

l, w = sm.symbols('l, w')
Ffl, Ffr, Frl, Frr = me.dynamicsymbols('F_{fl}, F_{fr}, F_{rl}, F_{rr}')
alphafl, alphafr = me.dynamicsymbols(r'\alpha_{fl}, \alpha_{fr}')
alpharl, alpharr = me.dynamicsymbols(r'\alpha_{rl}, \alpha_{rr}')
delta = me.dynamicsymbols('delta')

With the symbols defined, I use some auxiliary reference frames to establish the orientations with \(B\) being the car body, \(W\) being the steered front wheels, and the others to establish the direction of the force at each wheel.

B = me.ReferenceFrame('B')
W = me.ReferenceFrame('W')
FR = me.ReferenceFrame('F_R')
FL = me.ReferenceFrame('F_L')
RR = me.ReferenceFrame('R_R')
RL = me.ReferenceFrame('R_L')

W.orient_axis(B, delta, B.z)
FR.orient_axis(W, alphafr, W.z)
FL.orient_axis(W, alphafl, W.z)
RR.orient_axis(B, alpharr, B.z)
RL.orient_axis(B, alpharl, B.z)

The resultant of the forces expressed in the \(B\) frame is then:

R = Ffl*FL.x + Ffr*FR.x + Frl*RL.x + Frr*RR.x
R.express(B).simplify()

\[\displaystyle (F_{fl} \cos{\left(\alpha_{fl} + \delta \right)} + F_{fr} \cos{\left(\alpha_{fr} + \delta \right)} + F_{rl} \cos{\left(\alpha_{rl} \right)} + F_{rr} \cos{\left(\alpha_{rr} \right)})\hat{b}_x + (F_{fl} \sin{\left(\alpha_{fl} + \delta \right)} + F_{fr} \sin{\left(\alpha_{fr} + \delta \right)} + F_{rl} \sin{\left(\alpha_{rl} \right)} + F_{rr} \sin{\left(\alpha_{rr} \right)})\hat{b}_y\]

This resultant is bound to a line of action through \(B_o\). The associated couple is then calculated as the total moment about \(B_o\):

T = (me.cross(l/2*B.x - w/2*B.y, Ffl*FL.x) +
     me.cross(l/2*B.x + w/2*B.y, Ffr*FR.x) +
     me.cross(-l/2*B.x - w/2*B.y, Frl*RL.x) +
     me.cross(-l/2*B.x + w/2*B.y, Frr*RR.x))
T = T.express(B).simplify()
T

\[\displaystyle (\frac{\left(l \sin{\left(\alpha_{fl} + \delta \right)} + w \cos{\left(\alpha_{fl} + \delta \right)}\right) F_{fl}}{2} + \frac{\left(l \sin{\left(\alpha_{fr} + \delta \right)} - w \cos{\left(\alpha_{fr} + \delta \right)}\right) F_{fr}}{2} - \frac{\left(l \sin{\left(\alpha_{rl} \right)} - w \cos{\left(\alpha_{rl} \right)}\right) F_{rl}}{2} - \frac{\left(l \sin{\left(\alpha_{rr} \right)} + w \cos{\left(\alpha_{rr} \right)}\right) F_{rr}}{2})\hat{b}_z\]

Since we can always describe the forces acting on a rigid body as a resultant force and an associate torque of a couple, we will often take advantage of this simpler form for constructing models.

Specifying Forces and Torques

Forces are bound vectors, so we have to define their lines of action. This is best done by specifying the points on which each force acts, thus we will always use a vector and a point to fully describe the force. Methods and functions in SymPy Mechanics that make use of forces will typically require a tuple containing a point and a vector, for example the resultant force \(R^{S/B_o}\) acting on the mass center of the car would be specified like so:

Bo = me.Point('Bo')
force = (Bo, R)
force

(Bo, F_{fr}(t)*F_R.x + F_{fl}(t)*F_L.x + F_{rr}(t)*R_R.x + F_{rl}(t)*R_L.x)

Torques of a couple are free vectors (not bound to a line of action) but represent a couple acting on a rigid body, thus a reference frame associated with a rigid body and the vector representing the torque will be used to describe the torque in SymPy Mechanics. For example:

torque = (B, T)
torque

(B,
 ((l*sin(\alpha_{fl}(t) + delta(t)) + w*cos(\alpha_{fl}(t) + delta(t)))*F_{fl}(t)/2 + (l*sin(\alpha_{fr}(t) + delta(t)) - w*cos(\alpha_{fr}(t) + delta(t)))*F_{fr}(t)/2 - (l*sin(\alpha_{rl}(t)) - w*cos(\alpha_{rl}(t)))*F_{rl}(t)/2 - (l*sin(\alpha_{rr}(t)) + w*cos(\alpha_{rr}(t)))*F_{rr}(t)/2)*B.z)

We will often refer to forces and torques collectively as loads.

Note

The two cells above do not render the math nicely due to this SymPy bug: https://github.com/sympy/sympy/issues/24967.

Equal & Opposite

Both forces and torques applied to a multibody system must obey Newton’s Third Law, i.e. that forces and torques act equal and opposite. Take for example a torque from a motor that causes a pinned lever \(B\) to rotate relative to the ground \(N\) shown in Fig. 38. The motor torque can be modeled to occur between the stator and the rotor. We’ve arbitrarily selected the sign convention shown, i.e. a positive value of torque applies a positive torque to \(B\) and a negative torque to \(N\) if the torque is parallel to \(\hat{n}_z=\hat{b}_z\).

Fig. 38 A motor stator \(N\) fixed to ground with an arm fixed to the motor rotor \(B\) shown as one unit in a) and as separate bodies in b) with equal and opposite torque vectors applied to the pair of bodies representing the torque of a couple generated by the motor.

The motor torque can be specified as a time varying vector:

T, q = me.dynamicsymbols('T, q')

N = me.ReferenceFrame('N')
B = me.ReferenceFrame('B')

Tm = T*N.z

Then the equal and opposite torques are captured by these two tuples:

(B, Tm), (N, -Tm)

((B, T(t)*N.z), (N, - T(t)*N.z))

with equal and opposite torques applied to each body.

Warning

The sign conventions are really just a convention. It is also valid to choose (B, -Tm), (N, Tm) or even (B, Tm), (N, Tm) and (B, -Tm), (B, -Tm). But it is useful to choose a sign convention such that when the signs of angular velocity and torque are the same it corresponds to power into the system. So, for example, B.orient_axis(N, q, N.z) corresponds to (T*N.z, B) and power in. The key thing is that you know what your convention is so that you can interpret numerical results and signs correctly.

Contributing and Noncontributing Forces

Contributing forces are those that can do work on the multibody system. Work of a force \(\bar{F}\) acting over path \(S\) is defined as the following line integral:

(163)\[W = \int_S \bar{F} \cdot d\bar{s}\]

where \(d\bar{s}\) is the differential vector tangent to the path at the point the force is applied.

For example, the gravitational force acting on a particle moving through a unidirectional constant gravitational field (i.e. where the gravitational force is equal in magnitude, doesn’t change, and always the same direction) does work on the system unless the particle moves perpendicular to the field.

Noncontributing forces never do work on the system. For example, when a force acts between two points that have no relative motion, no work is done. Examples of noncontributing forces:

1. contact forces on particles across smooth (frictionless) surfaces of rigid bodies

2. any internal contact and body (distance) forces between any two points in a rigid body

3. contact forces between bodies rolling without slipping on each other

In the next chapter, we will see how the use of generalized coordinates relieve us from having to specify any noncontributing forces.

Gravity

We will often be interested in a multibody system’s motion when it is subject to gravitational forces. The simplest case is a constant unidirectional gravitational field, which is an appropriate model for objects moving on and near the Earth’s surface. The gravitational forces can be applied solely to the mass centers of each rigid body as a resultant force. The gravitational torque on the bodies is zero because the force is equal in magnitude for each particle in the body. See [Kane1985] pg. 110 for the more general model of Newton’s Law of Universal Gravitation where this is not the case. Studies of spacecraft dynamics often require considering both gravitational forces and moments.

In SymPy Mechanics, a gravitational force acting on a particle of mass \(m\) with acceleration due to gravity being \(g\) in the \(-\hat{n}_y\) direction would take this form:

m, g = sm.symbols('m, g')
Fg = -m*g*N.y
Fg

\[\displaystyle - g m\hat{n}_y\]

Springs & Dampers

Idealized springs and dampers are useful models of elements that have distance and velocity dependent forces and torques. A spring with free length \(q_0\) and where \(q_1,q_2\) locate the ends of the spring along a line parallel to \(\hat{n}_x\) is shown in Fig. 39.

If we displace \(P\) in the positive \(\hat{n}_x\) direction the spring will apply a force in the negative \(\hat{n}_x\) direction on point \(P\). So we chose a sign convention that the force on \(P\) from the spring is opposite the direction of the displacement.

Fig. 39 Diagram of a spring with a sign convention that tension is positive. \(P\) is shown separated from the end of the spring to show the equal and opposite forces.

If the spring is linear with stiffness \(k\) the spring force vector is then:

q0, k = sm.symbols('q0, k')
q1, q2 = me.dynamicsymbols('q1, q2')

displacement = q2 - q1 - q0
displacement

\[\displaystyle - q_{0} - q_{1} + q_{2}\]

Here a positive displacement represents the spring in tension and a negative displacement is compression.

Fs = -k*displacement*N.x
Fs

\[\displaystyle - k \left(- q_{0} - q_{1} + q_{2}\right)\hat{n}_x\]

Dampers are often used in parallel or series with springs to provide an energy dissipation via viscous-like friction. Springs combined with dampers allow for classical second order under-, over-, and critically-damped motion. A linear viscous damper with damping coefficient \(c\) that resists motion can be defined like so:

c = sm.symbols('c')
t = me.dynamicsymbols._t

Fc = -c*displacement.diff(t)*N.x
Fc

\[\displaystyle - c \left(- \dot{q}_{1} + \dot{q}_{2}\right)\hat{n}_x\]

Friction

Coulomb’s Law of Friction provides simple model of dry friction between two objects. When the two objects are in motion with respect to each other, there is a constant magnitude force that resists the motion. The force is independent of contact area and is proportional to the normal force between the objects. Coulomb’s kinetic friction model takes the scalar form:

(164)\[\begin{split}F_f = \begin{cases} \mu_k F_n & v < 0 \\ 0 & v = 0 \\ -\mu_k F_n & v > 0 \end{cases}\end{split}\]

where \(F_N\) is the normal force between the two objects, \(v\) is the relative speed between the two objects, and \(\mu_k\) is the coefficient of kinetic friction. At \(v=0\) kinetic friction is zero, but two objects in contact with a normal force can resist moving through static friction. When \(v=0\) any force perpendicular to the normal force can be generated up to a magnitude of \(F_f=\mu_s F_n\) where \(\mu_s\) is the coefficient of static friction and \(\mu_s > \mu_k\). Eq. (164) leaves this static case ambiguous, but it can be extended to:

(165)\[\begin{split}F_f = \begin{cases} \mu_k F_n & v < 0 \\ \left[-\mu_s F_n, \mu_s F_n\right] & v = 0 \\ -\mu_k F_n & v > 0 \end{cases}\end{split}\]

SymPy’s Piecewise is one way to create a symbolic representation of kinetic friction:

mu, m, g = sm.symbols('mu, m, g')

Fn = m*g

displacement = q2 - q1

Ff = sm.Piecewise((mu*Fn, displacement.diff(t) < 0),
                  (-mu*Fn, displacement.diff(t) > 0),
                  (0, True))*N.x
Ff

\[\begin{split}\displaystyle \begin{cases} g m \mu & \text{for}\: \dot{q}_{1} - \dot{q}_{2} > 0 \\- g m \mu & \text{for}\: \dot{q}_{1} - \dot{q}_{2} < 0 \\0 & \text{otherwise} \end{cases}\hat{n}_x\end{split}\]

The signum function (sign) can also be used in a similar and simpler form:

Ff = -mu*Fn*sm.sign(displacement.diff(t))*N.x
Ff

\[\displaystyle - g m \mu \operatorname{sign}{\left(- \dot{q}_{1} + \dot{q}_{2} \right)}\hat{n}_x\]

Eq. (165) is a sufficient model for many use cases, but it does not necessarily capture all observed effects. Fig. 40 shows a modification of Coulomb model that includes the Stribeck effect and viscous friction. Flores et. al have a nice summary of several other friction models that could be used [Flores2023].

Fig. 40 Extensions to the (a) Coulomb Dry Friction model: (b) Stribeck effect and (c) Stribeck and viscous effects. Taken from [Flores2023] (Creative Commons BY-NC-ND 4.0).

Aerodynamic Drag

Aerodynamic drag of a blunt body at low Reynolds numbers is dominated by the frontal area drag and the magnitude of this drag force can be modeled with the following equation:

(166)\[F_d = \frac{1}{2}\rho C_dAv^2\]

where \(\rho\) is the density of the air, \(C_d\) is the drag coefficient, \(A\) is the frontal area, and \(v\) is the air speed relative to the body.

If a body is moving in still air at an arbitrary velocity and point \(P\) is the aerodynamic center of the body then the aerodynamic drag force vector that opposes the motion can be found with such an equation:

A, Cd, rho = sm.symbols('A, C_d, rho')
ux, uy, uz = me.dynamicsymbols('u_x, u_y, u_z', real=True)

N_v_P = ux*N.x + uy*N.y + uz*N.z

Fd = -N_v_P.normalize()*Cd*A*rho/2*N_v_P.dot(N_v_P)
Fd

\[\displaystyle - \frac{A C_{d} \rho \sqrt{u_{x}^{2} + u_{y}^{2} + u_{z}^{2}} u_{x}}{2}\hat{n}_x - \frac{A C_{d} \rho \sqrt{u_{x}^{2} + u_{y}^{2} + u_{z}^{2}} u_{y}}{2}\hat{n}_y - \frac{A C_{d} \rho \sqrt{u_{x}^{2} + u_{y}^{2} + u_{z}^{2}} u_{z}}{2}\hat{n}_z\]

If the motion is only along the \(\hat{n}_x\) direction, for example, the equation for the drag force vector reduces to:

Fd.xreplace({uy: 0, uz:0})

\[\displaystyle - \frac{A C_{d} \rho u_{x} \left|{u_{x}}\right|}{2}\hat{n}_x\]

Managing the correct direction of the force, so that it opposes motion and is applied at the aerodynamic center, is important. The drag coefficient and frontal area can also change dynamically depending on the shape of the object and the direction the air is flowing over it.

Collision

If two points, a point and a surface, or two surfaces collide the impact behavior depends on the material properties, mass, and kinematics of the colliding bodies. There are two general approaches to modeling collision. The first is the Newtonian method in which you consider the momentum change, impulse, before and after collision. For a particle impacting a surface, this takes the basic form:

(167)\[m v^{+} = -e m v^{-}\]

where \(m\) is the particle’s mass, \(v^{-}\) is the speed before impact, \(v^{+}\) is the speed after impact, and \(e\) is the coefficient of restitution. The momentum after impact will be opposite and equal to the momentum before impact for a purely elastic collision \(e=1\) and the magnitude of the momentum will be less if the collision is inelastic \(0<e<1\). This approach can be extended to a multibody system; see [Flores2023] for an introduction to this approach.

The Newtonian model does not consider the explicit behavior of the force that generates the impulse at collision. Here we will take an alternative approach by modeling the force explicitly. Such contact force models can provide more accurate results, at the cost of longer computation times. Most impact force models build upon Hunt and Crossley’s seminal model [Hunt1975] which is based on Hertzian contact theory. Hunt and Crossley model the impact as a nonlinear function of penetration depth and its rate. The force is made up of a nonlinear stiffness and a damping term that take this form:

(168)\[f_c = k z^n + cz^n \dot{z}\]

\(k\) is the nonlinear contact stiffness, \(n\) is the stiffness exponent, \(z\) the contact penetration, \(\dot{z}\) is the penetration velocity, and \(c\) is the hysteresis damping factor. The damping scales with the penetration depth. \(k\) and \(c\) can be determined from the material properties and the shape of the colliding objects and can be related to the coefficient of restitution. \(n\) is \(3/2\) based on the Hertzian contact theory.

Fig. 41 Particle \(P\) colliding with a surface.

For example, if modeling a particle \(P\) that impacts a surface normal to \(\hat{n}_z\) that contains point \(O\) the penetration \(z_p\) of the particle into the surface (if positive \(z\) is out and negative \(z\) is inside the surface) can be described with:

(169)\[z_p = \frac{| \bar{r}^{P/O} \cdot \hat{n}_z | - \bar{r}^{P/O} \cdot \hat{n}_z}{2}\]

This difference between the absolute value and the value itself is equivalent to this piecewise function:

(170)\[\begin{split}z_p = \begin{cases} 0 & \bar{r}^{P/O} \cdot \hat{n}_z > 0 \\ \bar{r}^{P/O} \cdot \hat{n}_z & \bar{r}^{P/O} \cdot \hat{n}_z \leq 0 \end{cases}\end{split}\]

In SymPy, this can be defined like so:

x, y, z, zd = me.dynamicsymbols('x, y, z, \dot{z}', real=True)

r_O_P = x*N.x + y*N.y + z*N.z

zh = r_O_P.dot(N.z)

zp = (sm.Abs(zh) - zh)/2
zp

\[\displaystyle - \frac{z}{2} + \frac{\left|{z}\right|}{2}\]

The force can now be formulated according to (168):

k, c = sm.symbols('k, c')

Fz = (k*zp**(sm.S(3)/2) + c*zp**(sm.S(3)/2)*zd)*N.z
Fz

\[\displaystyle (c \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}} \dot{z} + k \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}})\hat{n}_z\]

We can check whether the force is correct for positive and negative \(z\):

Fz.xreplace({z: sm.Symbol('z', positive=True)})

\[\displaystyle 0\]

Fz.xreplace({z: sm.Symbol('z', negative=True)})

\[\displaystyle (c \left(- z\right)^{\frac{3}{2}} \dot{z} + k \left(- z\right)^{\frac{3}{2}})\hat{n}_z\]

More on the Hunt-Crossley model and alterations on the model are summarized in [Flores2023].

The impact force model is often combined with a friction model to generate a friction force for impacts that are not perfectly normal to the contacting surfaces. For example, Coulomb friction force can slow the particle’s sliding on the surface if we know the tangential velocity components \(v_x\) and \(v_y\) at the contact location. This lets us write to tangential friction force components:

mu = sm.symbols('mu')

vx = r_O_P.dot(N.x).diff(t)
vy = r_O_P.dot(N.y).diff(t)

Fx = -sm.Abs(vx)/vx*mu*Fz.dot(N.z)*N.x
Fx

\[\displaystyle - \frac{\mu \left(c \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}} \dot{z} + k \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}}\right) \left|{\dot{x}}\right|}{\dot{x}}\hat{n}_x\]

Fy = -sm.Abs(vy)/vy*mu*Fz.dot(N.z)*N.y
Fy

\[\displaystyle - \frac{\mu \left(c \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}} \dot{z} + k \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}}\right) \left|{\dot{y}}\right|}{\dot{y}}\hat{n}_y\]

These measure numbers for the force vector then evaluate to zero when there is no penetration \(z_p\) and evaluates to a spring and damper and Coulomb friction when there is. For example, using so numerical values to set the penetration:

vz = me.dynamicsymbols('v_z', negative=True)

repl = {zd: vz, z: sm.Symbol('z', positive=True)}

Fx.xreplace(repl), Fy.xreplace(repl), Fz.xreplace(repl)

\[\displaystyle \left( 0, \ 0, \ 0\right)\]

vz = me.dynamicsymbols('v_z', negative=True)

repl = {zd: vz, z: sm.Symbol('z', negative=True)}

Fx.xreplace(repl), Fy.xreplace(repl), Fz.xreplace(repl)

\[\displaystyle \left( - \frac{\mu \left(c \left(- z\right)^{\frac{3}{2}} v_{z} + k \left(- z\right)^{\frac{3}{2}}\right) \left|{\dot{x}}\right|}{\dot{x}}\hat{n}_x, \ - \frac{\mu \left(c \left(- z\right)^{\frac{3}{2}} v_{z} + k \left(- z\right)^{\frac{3}{2}}\right) \left|{\dot{y}}\right|}{\dot{y}}\hat{n}_y, \ (c \left(- z\right)^{\frac{3}{2}} v_{z} + k \left(- z\right)^{\frac{3}{2}})\hat{n}_z\right)\]

Finally, the total force on the particle contacting the surface can be fully described:

Fx + Fy + Fz

\[\displaystyle - \frac{\mu \left(c \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}} \dot{z} + k \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}}\right) \left|{\dot{x}}\right|}{\dot{x}}\hat{n}_x - \frac{\mu \left(c \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}} \dot{z} + k \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}}\right) \left|{\dot{y}}\right|}{\dot{y}}\hat{n}_y + (c \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}} \dot{z} + k \left(- \frac{z}{2} + \frac{\left|{z}\right|}{2}\right)^{\frac{3}{2}})\hat{n}_z\]