Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Function to return Total Energy in physics/mechanics #16202

Closed
wants to merge 13 commits into from
Closed
35 changes: 35 additions & 0 deletions sympy/physics/mechanics/particle.py
Expand Up @@ -223,3 +223,38 @@ def potential_energy(self, scalar):
"""

self._pe = sympify(scalar)

def total_energy(self, frame):
"""The total energy of the Particle.

The total energy, E, of a particle, P, is given by

'E = K.E + P.E.'

where P.E. is the potential energy of the particle P, and
K.E. is kinetic energy of the particle P, in the supplied ReferenceFrame.

Parameters
==========

frame : ReferenceFrame
The Particle's velocity is typically defined with respect to
an inertial frame but any relevant frame in which the velocity is
known can be supplied.

Examples
========

>>> from sympy.physics.mechanics import Particle, Point, ReferenceFrame
>>> from sympy import symbols
>>> m, v, r, g, h = symbols('m v r g h')
>>> N = ReferenceFrame('N')
>>> O = Point('O')
>>> P = Particle('P', O, m)
>>> P.point.set_vel(N, v * N.y)
>>> P.potential_energy = m * g * h
>>> P.total_energy(N)
g*h*m + m*v**2/2

"""
return self.kinetic_energy(frame) + self.potential_energy
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should raise an error if potential_energy is not defined.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At present, potential_energy for a body for which potential_energy is not explicitly defined is taken as zero. I think this makes sense. What do you think about it?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It is a bit misleading. If I type "total_energy" I'm expecting to get KE + PE, but if PE is zero by default I'm not getting any PE contribution. I feel like it would be best to inform the user that they aren't getting what they ask for. Maybe a warning message that tells the user that PE is zero and not set is sufficient.

128 changes: 121 additions & 7 deletions sympy/physics/mechanics/rigidbody.py
Expand Up @@ -199,6 +199,89 @@ def angular_momentum(self, point, frame):

return I.dot(w) + r.cross(m * v)

def rotational_kinetic_energy(self, frame):
"""Rotational kinetic energy of the rigid body

The rotational kinetic energy, T, of a rigid body, B, is given by

'T = 1/2 (I omega^2)'

where I is the central inertia dyadic of rigid body B,
omega is the angular velocity of the body's
mass center in the supplied ReferenceFrame.

Parameters
==========

frame : ReferenceFrame
The RigidBody's angular velocity and the velocity of it's mass
center are typically defined with respect to an inertial frame but
any relevant frame in which the velocities are known can be supplied.

Examples
========

>>> from sympy.physics.mechanics import Point, ReferenceFrame, outer
>>> from sympy.physics.mechanics import RigidBody
>>> from sympy import symbols
>>> M, v, r, omega = symbols('M v r omega')
>>> N = ReferenceFrame('N')
>>> b = ReferenceFrame('b')
>>> b.set_ang_vel(N, omega * b.x)
>>> P = Point('P')
>>> P.set_vel(N, v * N.x)
>>> I = outer (b.x, b.x)
>>> inertia_tuple = (I, P)
>>> B = RigidBody('B', P, b, M, inertia_tuple)
>>> B.rotational_kinetic_energy(N)
omega**2/2

"""

return (self.frame.ang_vel_in(frame) & (self.central_inertia &
self.frame.ang_vel_in(frame)) / sympify(2))
ritesh99rakesh marked this conversation as resolved.
Show resolved Hide resolved

def translational_kinetic_energy(self, frame):
"""Translational kinetic energy of the rigid body

The translational kinetic energy, T, of a rigid body, B, is given by

'T = 1/2 (m v^2)'

where m is the mass of rigid body B, v is the velocity of the body's
mass center in the supplied ReferenceFrame.

Parameters
==========

frame : ReferenceFrame
The RigidBody's angular velocity and the velocity of it's mass
center are typically defined with respect to an inertial frame but
any relevant frame in which the velocities are known can be supplied.

Examples
========

>>> from sympy.physics.mechanics import Point, ReferenceFrame, outer
>>> from sympy.physics.mechanics import RigidBody
>>> from sympy import symbols
>>> M, v, r, omega = symbols('M v r omega')
>>> N = ReferenceFrame('N')
>>> b = ReferenceFrame('b')
>>> b.set_ang_vel(N, omega * b.x)
>>> P = Point('P')
>>> P.set_vel(N, v * N.x)
>>> I = outer (b.x, b.x)
>>> inertia_tuple = (I, P)
>>> B = RigidBody('B', P, b, M, inertia_tuple)
>>> B.translational_kinetic_energy(N)
M*v**2/2

"""

return (self.mass * (self.masscenter.vel(frame) &
self.masscenter.vel(frame)) / sympify(2))
ritesh99rakesh marked this conversation as resolved.
Show resolved Hide resolved

def kinetic_energy(self, frame):
"""Kinetic energy of the rigid body

Expand Down Expand Up @@ -238,13 +321,7 @@ def kinetic_energy(self, frame):

"""

rotational_KE = (self.frame.ang_vel_in(frame) & (self.central_inertia &
self.frame.ang_vel_in(frame)) / sympify(2))

translational_KE = (self.mass * (self.masscenter.vel(frame) &
self.masscenter.vel(frame)) / sympify(2))

return rotational_KE + translational_KE
return self.rotational_kinetic_energy(frame) + self.translational_kinetic_energy(frame)

@property
def potential_energy(self):
Expand Down Expand Up @@ -296,3 +373,40 @@ def potential_energy(self, scalar):
"""

self._pe = sympify(scalar)

def total_energy(self, frame):
"""The total energy of the RigidBody.

The total energy E, of rigidbody B, is given by
E = K.E. + P.E.

where P.E. is the potential energy of the rigid body B and
K.E. is kinetic energy of the rigid body B in the supplied ReferenceFrame.

Parameters
==========

frame : ReferenceFrame
The RigidBody's angular velocity and the velocity of it's mass
center are typically defined with respect to an inertial frame but
any relevant frame in which the velocities are known can be supplied.

Examples
========

>>> from sympy.physics.mechanics import RigidBody, Point, outer, ReferenceFrame
>>> from sympy import symbols
>>> M, v, r, omega, g, h = symbols('M v r omega g h')
>>> N = ReferenceFrame('N')
>>> b = ReferenceFrame('b')
>>> b.set_ang_vel(N, omega * b.x)
>>> P = Point('P')
>>> P.set_vel(N, v * N.x)
>>> I = outer (b.x, b.x)
>>> Inertia_tuple = (I, P)
>>> B = RigidBody('B', P, b, M, Inertia_tuple)
>>> B.potential_energy = M * g * h
>>> B.total_energy(N)
M*g*h + M*v**2/2 + omega**2/2
"""
return self.kinetic_energy(frame) + self.potential_energy
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This should raise an error if potential energy is not defined.

6 changes: 3 additions & 3 deletions sympy/physics/mechanics/tests/test_particle.py
Expand Up @@ -33,7 +33,7 @@ def test_particle():
assert p.angular_momentum(O, N) == m2 * r * (v3 * N.x - v1 * N.z)
p.potential_energy = m * g * h
assert p.potential_energy == m * g * h
# TODO make the result not be system-dependent
assert p.total_energy(
N) == m*g*h + m2 * v1**2 / 2 + m2 * v2**2 / 2 + m2 * v3**2 / 2
assert p.kinetic_energy(
N) in [m2*(v1**2 + v2**2 + v3**2)/2,
m2 * v1**2 / 2 + m2 * v2**2 / 2 + m2 * v3**2 / 2]
N) == m2 * v1**2 / 2 + m2 * v2**2 / 2 + m2 * v3**2 / 2
ritesh99rakesh marked this conversation as resolved.
Show resolved Hide resolved
3 changes: 3 additions & 0 deletions sympy/physics/mechanics/tests/test_rigidbody.py
Expand Up @@ -53,7 +53,10 @@ def test_rigidbody2():
assert B.angular_momentum(O, N) == omega * b.x - M*v*r*b.z
B.potential_energy = M * g * h
assert B.potential_energy == M * g * h
assert expand(2 * B.translational_kinetic_energy(N)) == M * v**2
assert expand(2 * B.rotational_kinetic_energy(N)) == omega**2
assert expand(2 * B.kinetic_energy(N)) == omega**2 + M * v**2
assert B.total_energy(N) == M * g * h + M * v**2/2 + omega**2/2

def test_rigidbody3():
q1, q2, q3, q4 = dynamicsymbols('q1:5')
Expand Down