Add magnetic bearing feature

[tonylee2016]

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]

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]

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]

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]

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

[ross-bott]

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

• If the issue refers to a large task, break it in smaller issues that can be solved in
  less than 45 days;
• Label the issue as wontfix or wontfix for now and close it.

[mcarolalb]

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

[raphaeltimbo]

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

[mcarolalb]

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.

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:

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

[ross-bott]

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

• If the issue refers to a large task, break it in smaller issues that can be solved in
  less than 45 days;
• Label the issue as wontfix or wontfix for now and close it.