Repository: ross-rotordynamics/ross
Pull request #359 (2019-10-23 23:10:38)

Find eccentricity iteratively
Based on a given load, find eccentricity iteratively to match the equilibrium position.


Patch:
diff --git a/ross/fluid_flow/fluid_flow.py b/ross/fluid_flow/fluid_flow.py
index 34d38897..58b8229f 100644
--- a/ross/fluid_flow/fluid_flow.py
+++ b/ross/fluid_flow/fluid_flow.py
@@ -503,6 +503,35 @@ def fluid_flow_example():
     return my_pressure_matrix
 
 
+def fluid_flow_example2():
+    """This function returns a different instance of a simple fluid flow.
+    The purpose is to make available a simple model
+    so that doctest can be written using it.
+
+    Parameters
+    ----------
+
+    Returns
+    -------
+    An instance of a fluid flow object.
+
+    Examples
+    --------
+    >>> my_fluid_flow = fluid_flow_example2()
+    >>> my_fluid_flow.load
+    42.89
+    """
+    omega = 100. * 2 * np.pi / 60
+    radius_rotor = 0.1999996
+    radius_stator = 0.1999996 + 0.000194564
+    length = (1 / 10) * (2 * radius_stator)
+    viscosity = 0.015
+    density = 860
+    load = 42.89
+    my_pressure_matrix = FluidFlow(nz=8, ntheta=132, nradius=11, length=length, omega=omega,
+                                   p_in=0., p_out=0., radius_rotor=radius_rotor, radius_stator=radius_stator,
+                                   viscosity=viscosity, density=density, load=load)
+    return my_pressure_matrix
 
 
 
diff --git a/ross/fluid_flow/fluid_flow_coefficients.py b/ross/fluid_flow/fluid_flow_coefficients.py
index a138c53e..a76e1d5e 100644
--- a/ross/fluid_flow/fluid_flow_coefficients.py
+++ b/ross/fluid_flow/fluid_flow_coefficients.py
@@ -1,6 +1,7 @@
 import numpy as np
 from scipy import integrate
 from math import isnan
+from ross.fluid_flow.fluid_flow_geometry import calculate_rotor_load
 
 
 def calculate_analytical_stiffness_matrix(load, eccentricity_ratio, radial_clearance):
@@ -179,3 +180,80 @@ def calculate_oil_film_force(fluid_flow_object, force_type=None):
     force_y = radial_force * np.cos(fluid_flow_object.attitude_angle) \
         + tangential_force * np.sin(fluid_flow_object.attitude_angle)
     return radial_force, tangential_force, force_x, force_y
+
+
+def find_equilibrium_position(fluid_flow_object, print_along=True, relative_tolerance=1e-07,
+                              numerical_fit=False):
+    """This function returns an eccentricity value with calculated forces matching the load applied,
+    meaning an equilibrium position of the rotor.
+    Parameters
+    ----------
+    fluid_flow_object: A FluidFlow object.
+    print_along: bool, optional
+        If True, prints the iteration process.
+    relative_tolerance: float, optional
+    numerical_fit: bool, optional
+        If True, it makes sure the result matches numerically, changing attributes of the
+        FluidFlow object.
+    Returns
+    -------
+    float
+        Eccentricity of the equilibrium position.
+    Examples
+    --------
+    >>> from ross.fluid_flow.fluid_flow import fluid_flow_example2
+    >>> my_fluid_flow = fluid_flow_example2()
+    >>> find_equilibrium_position(my_fluid_flow, print_along=False) # doctest: +ELLIPSIS
+    0.00010000...
+    """
+    fluid_flow_object.calculate_pressure_matrix_numerical()
+    load = calculate_rotor_load(fluid_flow_object.radius_stator, fluid_flow_object.omega,
+                                fluid_flow_object.viscosity, fluid_flow_object.length,
+                                fluid_flow_object.radial_clearance, fluid_flow_object.eccentricity_ratio)
+    error = (load - fluid_flow_object.load)/fluid_flow_object.load
+    increment = 0.1
+    k = 0
+    eccentricity = fluid_flow_object.eccentricity
+    while np.abs(error) > relative_tolerance:
+        eccentricity = eccentricity + (eccentricity * increment)
+        eccentricity_ratio = eccentricity / fluid_flow_object.difference_between_radius
+        load = calculate_rotor_load(fluid_flow_object.radius_stator, fluid_flow_object.omega,
+                                    fluid_flow_object.viscosity, fluid_flow_object.length,
+                                    fluid_flow_object.radial_clearance, eccentricity_ratio)
+        new_error = (load - fluid_flow_object.load) / fluid_flow_object.load
+        if print_along:
+            print("Iteration " + str(k))
+            print("Eccentricity: " + str(eccentricity))
+            print("Load: " + str(load))
+            print("Error: " + str(new_error))
+        if error * new_error < 0:
+            if print_along:
+                print("Error changed sign. Changing sign of increment and reducing it.")
+            increment = -increment/10
+        elif abs(new_error) > abs(error):
+            if print_along:
+                print("Error was greater than previous one. Changing sign of increment and slightly "
+                      "reducing it.")
+            increment = -increment/5
+        error = new_error
+        k += 1
+        if print_along:
+            print()
+    return eccentricity
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+


Test patch:
diff --git a/ross/tests/test_fluid_flow.py b/ross/tests/test_fluid_flow.py
index db62ac99..44cd9a37 100644
--- a/ross/tests/test_fluid_flow.py
+++ b/ross/tests/test_fluid_flow.py
@@ -6,7 +6,7 @@ from bokeh.plotting import figure
 import matplotlib.pyplot as plt
 from numpy.testing import assert_allclose
 from ross.fluid_flow.fluid_flow_coefficients import calculate_analytical_damping_matrix,\
-    calculate_analytical_stiffness_matrix, calculate_oil_film_force
+    calculate_analytical_stiffness_matrix, calculate_oil_film_force, find_equilibrium_position
 from ross.fluid_flow.fluid_flow_graphics import plot_shape, plot_eccentricity, plot_pressure_theta,\
     plot_pressure_z, matplot_shape, matplot_eccentricity, matplot_pressure_theta, matplot_pressure_z
 
@@ -206,3 +206,10 @@ def test_matplotlib_plots():
     assert isinstance(matplot_pressure_z(bearing), ax_type)
     assert isinstance(matplot_shape(bearing), ax_type)
 
+
+def test_find_equilibrium_position():
+    bearing = fluid_flow_short_load()
+    eccentricity = find_equilibrium_position(bearing, print_along=False)
+    assert_allclose(eccentricity, (bearing.radius_stator - bearing.radius_rotor)*0.2663, rtol=0.0002)
+
+