petrobras / ross

ROSS is a library written in Python for rotordynamic analysis.
https://ross.readthedocs.io
Apache License 2.0
122 stars 95 forks source link

Add magnetic bearing feature #736

Open tonylee2016 opened 3 years ago

tonylee2016 commented 3 years ago

Glad to see a rotordynamics simulation package in python. I had a plan to do so when studied at Texas A&M.

I can contribute to improving the AMB module

  1. magnetic bearing model
  2. power electronics
  3. controller
rodrigomoliveira1 commented 3 years ago

Hello @tonylee2016! it's great to see your enthusiasm for contributing to ROSS.

Currently, we have a class to build a magnetic bearing. It's called MagneticBearingElement (see code below from ross/bearing_seal_element.py). You can take a step forward from what we have now. You're welcome to propose any improvement to the code.

@davidjgm is the one responsible for add this content to ROSS. I believe you both can help each other, if it's the case.

class MagneticBearingElement(BearingElement):
    """Magnetic bearing.
    This class creates a magnetic bearing element.
    Converts electromagnetic parameters and PID gains to stiffness and damping
    coefficients.
    Parameters
    ----------
    n : int
        The node in which the magnetic bearing will be located in the rotor.
    g0 : float
        Air gap in m^2.
    i0 : float
        Bias current in Ampere
    ag : float
        Pole area in m^2.
    nw : float or int
        Number of windings
    alpha : float or int
        Pole angle in radians.
    kp_pid : float or int
        Proportional gain of the PID controller.
    kd_pid : float or int
        Derivative gain of the PID controller.
    k_amp : float or int
        Gain of the amplifier model.
    k_sense : float or int
        Gain of the sensor model.
    tag : str, optional
        A tag to name the element
        Default is None.
    n_link : int, optional
        Node to which the bearing will connect. If None the bearing is
        connected to ground.
        Default is None.
    scale_factor : float, optional
        The scale factor is used to scale the bearing drawing.
        Default is 1.
    ----------
    See the following reference for the electromagnetic parameters g0, i0, ag, nw, alpha:
    Book: Magnetic Bearings. Theory, Design, and Application to Rotating Machinery
    Authors: Gerhard Schweitzer and Eric H. Maslen
    Page: 84-95
    Examples
    --------
    >>> n = 0
    >>> g0 = 1e-3
    >>> i0 = 1.0
    >>> ag = 1e-4
    >>> nw = 200
    >>> alpha = 0.392
    >>> kp_pid = 1.0
    >>> kd_pid = 1.0
    >>> k_amp = 1.0
    >>> k_sense = 1.0
    >>> tag = "magneticbearing"
    >>> mbearing = MagneticBearingElement(n=n, g0=g0, i0=i0, ag=ag, nw=nw,alpha=alpha,
    ...                                   kp_pid=kp_pid, kd_pid=kd_pid, k_amp=k_amp,
    ...                                   k_sense=k_sense)
    >>> mbearing.kxx
    [-4640.623377181318]
    """

    def __init__(
        self,
        n,
        g0,
        i0,
        ag,
        nw,
        alpha,
        kp_pid,
        kd_pid,
        k_amp,
        k_sense,
        tag=None,
        n_link=None,
        scale_factor=1,
    ):
        pL = [g0, i0, ag, nw, alpha, kp_pid, kd_pid, k_amp, k_sense]
        pA = [0, 0, 0, 0, 0, 0, 0, 0, 0]

        # Check if it is a number or a list with 2 items
        for i in range(9):
            if type(pL[i]) == float or int:
                pA[i] = np.array(pL[i])
            else:
                if type(pL[i]) == list:
                    if len(pL[i]) > 2:
                        raise ValueError(
                            "Parameters must be scalar or a list with 2 items"
                        )
                    else:
                        pA[i] = np.array(pL[i])
                else:
                    raise ValueError("Parameters must be scalar or a list with 2 items")

        # From: "Magnetic Bearings. Theory, Design, and Application to Rotating Machinery"
        # Authors: Gerhard Schweitzer and Eric H. Maslen
        # Page: 354
        ks = (
            -4.0
            * pA[1] ** 2.0
            * np.cos(pA[4])
            * 4.0
            * np.pi
            * 1e-7
            * pA[3] ** 2.0
            * pA[2]
            / (4.0 * pA[0] ** 3)
        )
        ki = (
            4.0
            * pA[1]
            * np.cos(pA[4])
            * 4.0
            * np.pi
            * 1e-7
            * pA[3] ** 2.0
            * pA[2]
            / (4.0 * pA[0] ** 2)
        )
        k = ki * pA[7] * pA[8] * (pA[5] + np.divide(ks, ki * pA[7] * pA[8]))
        c = ki * pA[7] * pA[5] * pA[8]
        # k = ki * k_amp*k_sense*(kp_pid+ np.divide(ks, ki*k_amp*k_sense))
        # c = ki*k_amp*kd_pid*k_sense

        # Get the parameters from k and c
        if np.isscalar(k):
            # If k is scalar, symmetry is assumed
            kxx = k
            kyy = k
        else:
            kxx = k[0]
            kyy = k[1]

        if np.isscalar(c):
            # If c is scalar, symmetry is assumed
            cxx = c
            cyy = c
        else:
            cxx = c[0]
            cyy = c[1]

        super().__init__(
            n=n,
            frequency=None,
            kxx=kxx,
            kxy=0.0,
            kyx=0.0,
            kyy=kyy,
            cxx=cxx,
            cxy=0.0,
            cyx=0.0,
            cyy=cyy,
            tag=tag,
            n_link=n_link,
            scale_factor=scale_factor,
        )
davidjgm commented 3 years ago

Hi @tonylee2016, Currently, the MagneticBearingElement class calculates the parameters of a PID controller based on the electromagnetic and geometric parameters. The reference book is: Magnetic Bearings. Theory, Design, and Application to Rotating Machinery, from Gerhard Schweitzer and Eric H. Maslen. Please, feel free to contribute to ROSS!

tonylee2016 commented 3 years ago

hi @davidjgm, have one question: I can not find the gain formulas used following on page 354 of Malsen's book, maybe a different version? can you let me know which chapter and section? https://github.com/ross-rotordynamics/ross/blob/84cc958b75be10bbf2473d2dd800762e0fb91730/ross/bearing_seal_element.py#L1508

davidjgm commented 3 years ago

Hi @tonylee2016, I miswrote the two references. The correct pages are: 69-80 and 343. I already commited the changes.

ross-bott commented 3 years ago

Hi there! I have marked this issue as stale because it has not had activity for 45 days. Consider the following options:

mcarolalb commented 1 year ago

Currently working on new ways to define the coeficients k and c of the magnetic bearing.

raphaeltimbo commented 4 months ago

As discussed with @ViniciusTxc3, new implementations for the magnetic bearing are being evaluated by @mcarolalb.

mcarolalb commented 2 months ago

Hi @raphaeltimbo!

In the updates regarding the implementation of a control system on ROSS, I conducted an initial test by creating a model of the magnetic bearing test rig that we have in our lab. The bearings are not shown in the figure because we will apply the magnetic force directly to the nodes where they are located, specifically nodes 12 and 43 of the model.

rotor

Then, I developed an initial code for the PID controller and the magnetic actuator.

    def PID1(self,Kp, Ki, Kd, setpoint, measurement, time):
        offset = 0
        e = setpoint - measurement
        P = Kp * e
        self.integral1 = self.integral1 + Ki * e * (time - self.time_prev1)
        D = Kd * (e - self.e_prev1) / (time - self.time_prev1)
        MV = offset + P + self.integral1 + D
        self.e_prev1 = e
        self.time_prev1 = time
        return MV
    def actuator(self,signal,displacement):
        const_k = (1/4) * (self.u0 * (self.n**2) * self.A)
        ki = (4 * const_k * self.i0) / (self.s0**2) * np.cos(self.alpha * (np.pi/180))
        ks = (-4 * const_k * (self.i0**2)) / (self.s0**3) * np.cos(self.alpha * (np.pi/180))
        magnetic_force = ki * signal + ks * displacement
        return magnetic_force

For the closed-loop controller system, I utilized the Newmark integrator and the add_to_RHS function to implement the magnetic force on the system.

    def controller(self, step, y, ydot):

        ext_force = np.zeros((rotor.ndof,))

        count = 0

        displacement = [y[4*self.node[0]], y[4*self.node[1]],y[4*self.node[0]+1], y[4*self.node[1]+1]]

        nos = [self.node[0], self.node[1], self.node[0], self.node[1]]

        for n, ii in zip(nos,displacement):            

            if count == 0:
                signal_pid1 = self.PID1(self.Kp, self.Ki, self.Kd, self.setpoint, ii, self.t[step])
                magnetic_force1 = self.atuador(signal_pid1,ii)
                ext_force[(4*n)+0] = magnetic_force1
            elif count == 1:
                signal_pid2 = self.PID2(self.Kp, self.Ki, self.Kd, self.setpoint, ii, self.t[step])
                magnetic_force2 = self.atuador(signal_pid2,ii)
                ext_force[(4*n)+0] = magnetic_force2
            elif count == 2:
                signal_pid3 = self.PID3(self.Kp, self.Ki, self.Kd, self.setpoint, ii, self.t[step])
                magnetic_force3 = self.atuador(signal_pid3,ii)
                ext_force[(4*n)+1] = magnetic_force3
            elif count == 3:
                signal_pid4 = self.PID4(self.Kp, self.Ki, self.Kd, self.setpoint, ii, self.t[step])
                magnetic_force4 = self.atuador(signal_pid4,ii)
                ext_force[(4*n)+1] = magnetic_force4

            count = count + 1

        return ext_force
        response = rotor.run_time_response(speed, F, t, integrator = "newmark", add_to_RHS = self.controller,progress_interval=1)

With this, I obtained the following results:

probe1

probe2

Now, I’m working on the final implementation for a better integration with ROSS.

ross-bott commented 3 weeks ago

Hi there! I have marked this issue as stale because it has not had activity for 45 days. Consider the following options: