diff --git a/CHANGELOG.md b/CHANGELOG.md
index 0bc6754f5..e682e8890 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -43,6 +43,7 @@ Attention: The newest changes should be on top -->
 
 ### Fixed
 
+- BUG: Zero Mass Flow Rate in Liquid Motors breaks Exhaust Velocity [#677](https://github.com/RocketPy-Team/RocketPy/pull/677)
 - DOC: Fix documentation dependencies [#651](https://github.com/RocketPy-Team/RocketPy/pull/651)
 - DOC: Fix documentation warnings [#645](https://github.com/RocketPy-Team/RocketPy/pull/645)
 - BUG: Rotational EOMs Not Relative To CDM [#674](https://github.com/RocketPy-Team/RocketPy/pull/674)
diff --git a/rocketpy/motors/liquid_motor.py b/rocketpy/motors/liquid_motor.py
index 09e3aadd5..9ec3d1130 100644
--- a/rocketpy/motors/liquid_motor.py
+++ b/rocketpy/motors/liquid_motor.py
@@ -266,16 +266,16 @@ def exhaust_velocity(self):
         """
         times, thrusts = self.thrust.source[:, 0], self.thrust.source[:, 1]
         mass_flow_rates = self.mass_flow_rate(times)
+        exhaust_velocity = np.zeros_like(mass_flow_rates)
 
         # Compute exhaust velocity only for non-zero mass flow rates
         valid_indices = mass_flow_rates != 0
-        valid_times = times[valid_indices]
-        valid_thrusts = thrusts[valid_indices]
-        valid_mass_flow_rates = mass_flow_rates[valid_indices]
 
-        ext_vel = -valid_thrusts / valid_mass_flow_rates
+        exhaust_velocity[valid_indices] = (
+            -thrusts[valid_indices] / mass_flow_rates[valid_indices]
+        )
 
-        return np.column_stack([valid_times, ext_vel])
+        return np.column_stack([times, exhaust_velocity])
 
     @funcify_method("Time (s)", "Propellant Mass (kg)")
     def propellant_mass(self):
diff --git a/rocketpy/motors/motor.py b/rocketpy/motors/motor.py
index 07d969e3d..eae7931b4 100644
--- a/rocketpy/motors/motor.py
+++ b/rocketpy/motors/motor.py
@@ -483,7 +483,8 @@ def total_mass_flow_rate(self):
         rate should not be greater than `total_mass_flow_rate`, otherwise the
         grains mass flow rate will be negative, losing physical meaning.
         """
-        return -1 * self.thrust / self.exhaust_velocity
+        average_exhaust_velocity = self.total_impulse / self.propellant_initial_mass
+        return self.thrust / -average_exhaust_velocity
 
     @property
     @abstractmethod