Diseño de un cohete de agua

In [56]:

    

from IPython.display import Image, SVG, Math, Latex, clear_output
import sys
import numpy as np

P0 = 1E5  # Pa
G = -9.81  # m/s
RHO_W = 998.0  # kg / m3
RHO_A = 1.0  # kg / m3
DT = 0.00001  # s
DT_BIG = 0.1  # s

# Constest related data
__L_WING__ = 0.025  # m
__Z_FS__ = 0.15  # m
__P__ = 7E5  # Pa
__ANGLE__ = np.radians(80.0)
__NU_PROP__ = 0.9


    

In [29]:

    

SVG(filename='./forces.svg')


    

    
Out[29]:

In [30]:

    

SVG(filename='./stable_unstable.svg')


    

    
Out[30]:

In [31]:

    

Image(filename='unstable_rocket.jpeg')


    

    
Out[31]:

In [32]:

    

SVG(filename='./scheme.svg')


    

    
Out[32]:

In [33]:

    

%pylab inline

# From http://www.mzbn.com/en/sales/information
R = 0.5 * 65E-3  # m
L = 223E-3  # m
# From http://www.embotelladoras.org/maquinas-llenadoras/llenadora-de-botella-de-agua-embotellada-1500-botellashora
r0_bottle = 0.5 * 25E-3  # m
r0_valve = r0_bottle / 3.0  # m

z = (0.0, 0.01, 0.02, 0.04, 0.06, 0.175, 0.195)
r = (r0_valve, r0_bottle, r0_bottle, 0.029, 0.036, 0.036, 0.01)
geom_z = np.asarray(z)
geom_r = np.asarray(r)

plot(geom_r, geom_z, 'k-')
plot(-geom_r, geom_z, 'k-')
plot((-geom_r[-1], geom_r[-1]), (geom_z[-1], geom_z[-1]), 'k-')
axes().set_aspect('equal', 'datalim')


    

    




Populating the interactive namespace from numpy and matplotlib






    

In [34]:

    

class areaElement(object):
    def __init__(self, z, r):
        """Initialize the rocket bottle geometry.
        
        Args:
            z: Points where the radius is provided
            r: Radius associated to each z points
        Raises:
            ValueError: if r and z have different length, or if the
            lengths are lower than 2.
        """
        if len(z) != len(r):
            self.__z = None
            self.__r = None
            raise ValueError('Unmatching dimensions for z and r')
        if len(z) < 2:
            self.__z = None
            self.__r = None
            raise ValueError('At least 2 points should be provided')            
        z = np.asarray(z)
        r = np.asarray(r)
        self.__z = z
        self.__r = r
        # Compute the volumes
        self.__V = np.pi * np.asarray(
            [np.trapz(r[0:i + 1]**2, x=z[0:i + 1]) for i in range(len(z))])
        # Compute the areas
        self.__A = np.pi * np.asarray(
            [np.trapz(r[0:i + 1], x=z[0:i + 1]) for i in range(len(z))])
        self.__S = np.pi * r**2
        # The centers of gravity
        self.__C = np.pi * np.asarray(
            [np.trapz(z[0:i + 1] * r[0:i + 1]**2, x=z[0:i + 1]) for i in range(len(z))])
        self.__C[1:] /= self.__V[1:]  # V[0] = 0.0

    @property
    def r(self):
        """List of radius R(z)"""
        return self.__r

    @r.setter
    def r(self, value):
        self.__r = value

    @property
    def z(self):
        """List of z values"""
        return self.__z

    @z.setter
    def z(self, value):
        self.__z = value

    def getR(self, z=None):
        """Get the radius of the bottle for a specific heigh.
        
        Args:
            z: Heigh where the radius should be known. If None, the upper bound
            of the bottle will be used.
        Return:
            The radius at z. z value will be clamped by the bottle bounds.
        """
        if z == None:
            return self.__r[-1]
        return np.interp(z, self.__z, self.__r)

    def getN(self, z=None):
        """Get the normal (lateral wall) of the bottle for a specific heigh.
        
        Args:
            z: Heigh where the normal should be known. If None, the upper bound
            of the bottle will be used.
        Return:
            The normal at z (nr, nz). z value will be clamped by the bottle bounds.
        """
        if z == None:
            return np.asarray((self.__nr[-1], self.__nz[-1]))
        nr = np.interp(z, self.__z, self.__nr)
        nz = np.interp(z, self.__z, self.__nz)
        return np.asarray((nr, nz))

    def getA(self, z=None):
        """Get the lateral area of the bottle below a especific z coordinate.
        
        Args:
            z: Heigh where the area should be known. If None, the upper bound
            of the bottle will be used.
        Return:
            The lateral area below z. z value will be clamped by the bottle bounds.
        """
        if z == None:
            return self.__A[-1]
        return np.interp(z, self.__z, self.__A)

    def getS(self, z=None):
        """Get the transversal area of the bottle for a specific heigh.
        
        Args:
            z: Heigh where the transversal area should be known. If None, the upper bound
            of the bottle will be used.
        Return:
            The lateral area below z. z value will be clamped by the bottle bounds.
        """
        if z == None:
            return self.__S[-1]
        return np.interp(z, self.__z, self.__S)

    def getV(self, z=None):
        """Get the Volume up to a specific heigh.
        
        Args:
            z: Upper bound of the volume. If None, the upper bound of the bottle
            will be used.
        Return:
            The volume under z. z value will be clamped to the bottle bounds.
        """
        if z == None:
            return self.__V[-1]
        return np.interp(z, self.__z, self.__V)

    def getC(self, z=None):
        """Get the Center of gravity of the volume getV().
        
        Args:
            z: Upper bound of the volume. If None, the upper bound of the bottle
            will be used.
        Return:
            The center of gravity of the volume under z.
        """
        if z == None:
            return self.__C[-1]
        return np.interp(z, self.__z, self.__C)

bottle_test = areaElement(geom_z, geom_r)
A = bottle_test.getA() + bottle_test.getS()
V = bottle_test.getV()
C = bottle_test.getC()
Math('A = {:.2g} \\, \\mathrm{{m^2}}; V = {:.2g} \\, \\mathrm{{m^3}}, C = {:.2g} \\, \\mathrm{{m}}'.format(
    A, V, C))


    

    
Out[34]:





$$A = 0.019 \, \mathrm{m^2}; V = 0.00062 \, \mathrm{m^3}, C = 0.11 \, \mathrm{m}$$

In [35]:

    

SVG(filename='./control_vol.svg')


    

    
Out[35]:

In [36]:

    

class structure(areaElement):
    def __init__(self, z, r):
        """Initialize the rocket bottle geometry.
        
        Args:
            z: Points where the radius is provided
            r: Radius associated to each z points
        Raises:
            ValueError: if r and z have different length, or if the
            lengths are lower than 2.
        """
        areaElement.__init__(self, z, r)

    def getF(self, p, z_FS, z=None):
        """Get the pressure integral up to a specific heigh.
        
        Args:
            p: Air pressure
            z_FS: z coordinate of the Free surface
            z: Heigh where the radius should be known. If None, the upper bound
            of the bottle will be used.
        Return:
            The vertical component of the pressure integral.
        """
        if z == None:
            z = self.z[-1]
        if z_FS <= 0.0 or z <= 0.0 or p == P0:
            return 0.0
        # Create the list of walls to integrate
        i = 0
        zz = [self.z[i]]
        rr = [self.r[i]]
        while self.z[i + 1] < z:
            i += 1
            zz.append(self.z[i])
            rr.append(self.r[i])
        zz.append(z)
        rr.append(self.getR(z))
        # Solve the taps integrals
        result = p * self.getS(z) - P0 * self.getS(0.0)
        # Integrate the lateral walls
        for i in range(1, len(zz)):
            z0 = zz[i - 1]
            z1 = zz[i]
            r0 = rr[i - 1]
            r1 = rr[i]
            zi = 0.5 * (z0 + z1)
            ri = 0.5 * (r0 + r1)
            dr = r1 - r0
            if zi > z_FS:
                press = p
            else:
                press = P0 + (p - P0) / z_FS * zi
            result -= 2.0 * np.pi * press * ri * dr
        return result

bottle_test = structure(geom_z, geom_r)
F = bottle_test.getF(p=7E5, z_FS=0.15)
Math(r'F = {:.2g} N'.format(F))


    

    
Out[36]:





$$F = 1.9e+03 N$$

In [37]:

    

class water(object):
    """Water volume inside the bottle.
    
    The water volume is deflated while the rocket is propelled, pushed by the
    air.
    """
    def __init__(self, bottle, z, rho=RHO_W):
        """Init the water volume.

        Args:
            bottle: Bottle where the water should be contained.
            z: Initial water level.
            rho: Water density
            
        """
        self.__bottle = bottle
        self.z = z
        self.rho = rho
        self.u = 0.0
        self.dudt = 0.0

    def update(self, p, dUdt, g=G, dt=DT):
        """Update the water free surface position.

        Args:
            p: Air pressure.
            dUdt: Rocket acceleration
            dt: Time step
        
        """
        z = self.z
        if z <= 0.0:
            return
        # Free surface acceleration
        rho = self.rho
        rz = self.__bottle.getR(z)
        r0 = self.__bottle.getR(0)
        V = self.getV()
        dudt = -rho * np.pi * (rz / r0)**2 * self.u**2 * (r0**2 - rz**2)
        dudt -= self.__bottle.getF(p=p, z_FS=z, z=z)
        dudt += rho * (g - dUdt) * V
        dudt /= rho * np.pi * rz**2 * z
        # Time integration
        self.u += dt * dudt
        self.z += dt * self.u
        self.dudt = dudt
        if self.z < 0.0:
            self.z = 0.0

    def getZ(self):
        """Compute the water level.

        Return:
            Level of water
        
        """
        return self.z

    def getV(self):
        """Compute the volume of water.

        Return:
            Volume of water
        
        """
        return self.__bottle.getV(self.z)

    def getM(self):
        """Compute the water mass.

        Return:
            Mass of water
        
        """
        return self.getV() * self.rho

    def getC(self):
        """Get the Center of gravity of the water mass.
        
        Return:
            The center of gravity of the water mass.
        """
        return self.__bottle.getC(self.z)


    

In [38]:

    

class air(object):
    """Air volume inside the bottle.
    
    The water volume is inflated while the rocket is propelled, pushing the
    water.
    """
    def __init__(self, bottle, z, p):
        """Init the air volume.

        Args:
            bottle: Bottle where the water should be contained.
            z: Initial water level.
            p: Initial air pressure
            
        """
        self.__bottle = bottle
        self.pV = p * (bottle.getV() - bottle.getV(z))
        self.p = p

    def update(self, z):
        """Update the air volume and pressure.

        Args:
            z: New free surface position.
        
        """
        if z <= 0.0:
            self.p = P0
            return
        V = self.__bottle.getV() - self.__bottle.getV(z)
        self.p = self.pV / V
        if self.p < P0:
            self.p = P0


    

In [39]:

    

class propeller(structure):
    def __init__(self, z, r, z_FS, p):
        """Initialize the rocket bottle geometry.
        
        Args:
            z: Points where the radius is provided
            r: Radius associated to each z points
            z_FS: Initial water level
            p: Initial pressure
        Raises:
            ValueError: if r and z have different length, or if the
            lengths are lower than 2.
        """
        structure.__init__(self, z, r)
        self.water = water(self, z_FS)
        self.air = air(self, z_FS, p)

    def update(self, dUdt, g=G, dt=DT):
        """Update the physical state.

        Args:
            dUdt: Acceleration of the rocket (forward pointing)
            dt: Time step.
        
        """
        # Update the water and air
        self.water.update(self.air.p, dUdt, g=g, dt=DT)
        self.air.update(self.water.z)

    def getFprop(self):
        """Get the propeller force.
        
        Return:
            The propeller force.
        """
        p = self.air.p
        z_FS = self.water.z
        return self.getF(p, z_FS)

z_FS = __Z_FS__
p = __P__
bottle_test = propeller(geom_z, geom_r, z_FS, p)
t = [0.0]
p = [p]
while bottle_test.water.z > 0.0:
    bottle_test.update(0.0, dt=DT)
    t.append(t[-1] + DT)
    p.append(bottle_test.air.p)
plot(t, p, 'k-')
grid()
xlabel(r'$t$ [s]')
ylabel(r'$p$ [Pa]')


    

    
Out[39]:





<matplotlib.text.Text at 0x7f1f07c2c128>






    

In [40]:

    

SVG(filename='./hulls.svg')


    

    
Out[40]:

In [41]:

    

Image(filename='torpedoe.jpeg')


    

    
Out[41]:

In [42]:

    

Image(filename='bullet.jpeg')


    

    
Out[42]:

In [43]:

    

Image(filename='missil.jpeg')


    

    
Out[43]:

In [44]:

    

Image(filename='missil_15_plot.png')


    

    
Out[44]:

In [45]:

    

class hull(structure):
    def __init__(self, z, r, alpha, CL, CD, Cz, rho=RHO_A):
        """Initialize the rocket bottle geometry.
        
        Args:
            z: Points where the radius is provided
            r: Radius associated to each z points
            alpha: List of angles from which the aerodynamic coefficients are provided (radians).
            CL: Lift coefficient for each angle.
            CD: Drag coefficient for each angle.
            Cz: aerodynamic center heigh coeffcient (addimensionalized with the cylinder heigh)
            rho: Fluid where the hull is moving
        Raises:
            ValueError: if r and z have different length, or if the
            lengths are lower than 2, or if alpha, CL, CD, Cz have different lengths,
            or if alpha, CL, CD, Cz are lower than 2.
        """
        structure.__init__(self, z, r)
        if len(alpha) != len(CL) != len(CD) != len(Cz):
            raise ValueError('Unmatching dimensions for alpha, CL, CD, Cz')
        if len(alpha) < 2:
            raise ValueError('At least two angles should be provided')

        self.__alpha = []
        self.__CL = []
        self.__CD = []
        self.__Cz = []
        for i in range(0, len(alpha)):
            self.__alpha.append(-alpha[-(i + 1)])
            self.__CL.append(-CL[-(i + 1)])
            self.__CD.append(CD[-(i + 1)])
            self.__Cz.append(Cz[-(i + 1)])
        for i in range(0, len(alpha)):
            self.__alpha.append(alpha[i])
            self.__CL.append(CL[i])
            self.__CD.append(CD[i])
            self.__Cz.append(Cz[i])
        self.__alpha = np.asarray(self.__alpha)
        self.__CL = np.asarray(self.__CL)
        self.__CD = np.asarray(self.__CD)
        self.__Cz = np.asarray(self.__Cz)
        self.__rho = rho

        # Compute the area
        self.__Area = self.getA() + self.getS() + self.getS(z[0])
        
        # Compute the cylinder heigh
        r_max = max(r)
        z_id = numpy.where(r == r_max)[0][-1]
        self.__z_cyl = z[z_id]
            
            
    def getL(self, alpha, U):
        """Compute the lift force.

        Args:
            alpha: Angle with the input flow (radians)
            U: Input flow velocity (equal to the rocket speed).
        Return:
            The lift force.
        """
        CL = np.interp(alpha, self.__alpha, self.__CL)
        return 0.5 * self.__rho * self.__Area * U**2 * CL

    def getD(self, alpha, U):
        """Compute the drag force.

        Args:
            alpha: Angle with the input flow (radians)
            U: Input flow velocity (equal to the rocket speed).
        Return:
            The drag force.
        """
        CD = np.interp(alpha, self.__alpha, self.__CD)
        return 0.5 * self.__rho * self.__Area * U**2 * CD

    def getZ(self, alpha, U):
        """Compute the lift and drag forces application point heigh.

        Args:
            alpha: Angle with the input flow (radians)
        Return:
            The application point heigh.
        """
        Cz = np.interp(alpha, self.__alpha, self.__Cz)
        return self.__z_cyl * Cz

z = (0.0, 0.175, 0.18444, 0.19444, 0.20421, 0.21255, 0.22042, 0.22723, 0.23306, 0.23724)
r = (0.036, 0.036, 0.03531, 0.03326, 0.02993, 0.02546, 0.02, 0.01378, 0.00702, 0.0)
hull_z = np.asarray(z)
hull_r = np.asarray(r)

plot(r, z, 'k-')
axes().set_aspect('equal', 'datalim')


    

In [46]:

    

alpha = np.asarray((np.radians(0.0), np.radians(5.0), np.radians(10.0), np.radians(15.0)))
CL = np.asarray((0.0, 0.016, 0.029, 0.035))
CD = np.asarray((0.023, 0.028, 0.033, 0.041))
Cz = np.asarray((0.68, 0.91, 0.98, 1.1))

hull_test = hull(hull_z, hull_r, alpha, CL, CD, Cz)
a = np.linspace(-15.0, 15.0, num=30)
l = []
d = []
z = []
for aa in a:
    l.append(hull_test.getL(np.radians(aa), 10.0))
    d.append(hull_test.getD(np.radians(aa), 10.0))
    z.append(hull_test.getZ(np.radians(aa), 10.0))
plot(a, l, 'r-', label='Lift')
plot(a, d, 'b-', label='Drag')
plot(a, z, 'k-', label='z')
legend(loc='best')
grid()
xlabel(r'$\alpha$ [deg]')
ylabel(r'$L, D, z$')


    

    
Out[46]:





<matplotlib.text.Text at 0x7f1f07c12438>






    

In [47]:

    

SVG(filename='./wing_scheme.svg')


    

    
Out[47]:

In [48]:

    

class wing():
    def __init__(self, z, r, alpha, CLo, CDo, rho=RHO_A):
        """Initialize the rocket bottle geometry.
        
        Args:
            z: Points where the radius is provided
            r: Radius associated to each z points
            alpha: List of angles from which the aerodynamic coefficients are provided (radians).
            CLo: Lift coefficient for each angle (without downwash effect).
            CDo: Drag coefficient for each angle (without downwash effect).
            rho: Fluid where the hull is moving
        Raises:
            ValueError: if r and z have different length, or if the
            lengths are lower than 2, or if alpha, CL, CD, Cz have different lengths,
            or if alpha, CL, CD, Cz are lower than 2.
        """
        if len(z) != len(r):
            self.__z = None
            self.__r = None
            raise ValueError('Unmatching dimensions for z and r')
        if len(z) < 2:
            self.__z = None
            self.__r = None
            raise ValueError('At least 2 points should be provided')
        z = np.asarray(z)
        r = np.asarray(r)
        self.__z = z
        self.__r = r

        if len(alpha) != len(CLo) != len(CDo):
            raise ValueError('Unmatching dimensions for alpha, CL, CD')
        if len(alpha) < 2:
            raise ValueError('At least two angles should be provided')
    
        self.__alpha = []
        self.__CLo = []
        self.__CDo = []
        for i in range(0, len(alpha)):
            self.__alpha.append(-alpha[-(i + 1)])
            self.__CLo.append(-CLo[-(i + 1)])
            self.__CDo.append(CDo[-(i + 1)])
        for i in range(0, len(alpha)):
            self.__alpha.append(alpha[i])
            self.__CLo.append(CLo[i])
            self.__CDo.append(CDo[i])
        self.__alpha = np.asarray(self.__alpha)
        self.__CLo = np.asarray(self.__CLo)
        self.__CDo = np.asarray(self.__CDo)
        self.__rho = rho

        # Compute the area
        self.__Area = np.trapz(r, x=z)

        # Compute the center of gravity
        self.__cor = np.trapz(r * z, x=z) / self.__Area

        # Compute the aspect ratio
        r_max = max(r)
        self.__Ar = r_max**2 / self.__Area

    @property
    def r(self):
        """List of radius R(z)"""
        return self.__r

    @r.setter
    def r(self, value):
        self.__r = value

    @property
    def z(self):
        """List of z values"""
        return self.__z

    @z.setter
    def z(self, value):
        self.__z = value
        
    def getL(self, alpha, U):
        """Compute the lift force.

        Args:
            alpha: Angle with the input flow (radians)
            U: Input flow velocity (equal to the rocket speed).
        Return:
            The lift force.
        """
        CLo = np.interp(alpha, self.__alpha, self.__CLo)
        CL = CLo / (1 + CLo / (np.pi * self.__Ar))
        return 0.5 * self.__rho * self.__Area * U**2 * CL

    def getD(self, alpha, U):
        """Compute the drag force.

        Args:
            alpha: Angle with the input flow (radians)
            U: Input flow velocity (equal to the rocket speed).
        Return:
            The drag force.
        """
        CDo = np.interp(alpha, self.__alpha, self.__CDo)
        CD = CDo + self.getL(alpha, U)**2 / (0.7 * np.pi * self.__Ar)
        return 0.5 * self.__rho * self.__Area * U**2 * CD

    def getZ(self, alpha, U):
        """Compute the lift and drag forces application point heigh.

        Args:
            alpha: Angle with the input flow (radians)
        Return:
            The application point heigh.
        """
        return self.__cor


z = (0.0, 0.006)
r = (__L_WING__, __L_WING__)
wing_z = np.asarray(z)
wing_r = np.asarray(r)

alpha = np.asarray((np.radians(0.0), np.radians(5.0), np.radians(10.0), np.radians(15.0)))
CLo = np.asarray((0.0, 0.5438, 0.9067, 0.4365))
CDo = np.asarray((0.0101, 0.0141, 0.0270, 0.1754))

wing_test = wing(wing_z, wing_r, alpha, CLo, CDo)
a = np.linspace(-15.0, 15.0, num=30)
l = []
d = []
z = []
for aa in a:
    l.append(wing_test.getL(np.radians(aa), 10.0))
    d.append(wing_test.getD(np.radians(aa), 10.0))
    z.append(wing_test.getZ(np.radians(aa), 10.0))
plot(a, l, 'r-', label='Lift')
plot(a, d, 'b-', label='Drag')
plot(a, z, 'k-', label='z')
legend(loc='best')
grid()
xlabel(r'$\alpha$ [deg]')
ylabel(r'$L, D, z$')


    

    
Out[48]:





<matplotlib.text.Text at 0x7f1f07a389b0>






    

In [49]:

    

class solid(object):
    def __init__(self, P, N, U=np.asarray((0.0, 0.0)), W=0.0):
        """Initialize the solid.
        
        Args:
            P: Initial Point
            N: Initial Direction
            U: Initial Velocity
            W: Initial Angular Velocity
        Raises:
            ValueError: if P, N, or U are not 2D entities.
        """
        if len(P) != len(N) != len(U) != 2:
            raise ValueError('P, N, and U must be 2D entities')
        self.__P = np.asarray(P)
        self.__N = np.asarray(N)
        self.__U = np.asarray(U)
        self.__W = W

    @property
    def P(self):
        """Position"""
        return self.__P

    @P.setter
    def P(self, value):
        self.__P = value

    @property
    def N(self):
        """Direction N"""
        return self.__N

    @N.setter
    def N(self, value):
        self.__N = value

    @property
    def U(self):
        """Velocity"""
        return self.__U

    @U.setter
    def U(self, value):
        self.__U = value

    @property
    def W(self):
        """Angular Velocity"""
        return self.__W

    @W.setter
    def W(self, value):
        self.__W = value

    def update(self, m, I, z_cog, F, z_F, dt=DT):
        """Time integration.
        
        Args:
            m: Rocket mass.
            I: Rocket inertia.
            z_cog: heigh of the center of gravity in local coordinates.
            F: List of forces applied to the rocket (in local coordinates).
            z_F: heigh of the centers of application of the forces in local coordinates.
        """
        # Compute the global force and moment
        Fl = np.asarray((0.0, 0.0))
        Ml = 0.0
        for i, f in enumerate(F):
            if len(f) != 2:
                raise ValueError('The forces should be 2D entities')
            f = np.asarray(f)
            Fl += f
            z_f = z_F[i] - z_cog
            Ml -= f[0] * z_f

        # Change the force to the global coordinates
        T = np.asarray((self.N[1], -self.N[0]))
        Fg = Fl[1] * self.N + Fl[0] * T
        Mg = Ml

        # Compute the acelerations
        dUdt = Fg / m
        dWdt = Mg / I

        # Integrate in time
        self.U += dt * dUdt
        self.P += dt * self.U
        self.W += dt * dWdt
        dTheta = dt * self.W
        self.N = np.matrix([[np.cos(dTheta), -np.sin(dTheta)],
                            [np.sin(dTheta), np.cos(dTheta)]]).dot(self.N)
        self.N = np.asarray(self.N)[0]
        return dUdt


    

In [54]:

    

class rocket(solid):
    def __init__(self, angle, m, I, z_cog, p, h, w):
        """Initialize the rocket.
        
        Args:
            angle: Launching angle (radians). 0 for the x axis.
            m: Rocket mass (without water)
            I: Rocket inertia (without water)
            z_cog: Rocket center of gravity heigh (without water)
            p: Propeller
            h: hull
            w: List of wings
        """
        # Solid initialization
        P = np.matrix([[np.cos(angle), -np.sin(angle)],
                       [np.sin(angle), np.cos(angle)]]).dot(np.asarray((z_cog, 0.0)))
        P = np.asarray(P)[0]
        N = np.matrix([[np.cos(angle), -np.sin(angle)],
                       [np.sin(angle), np.cos(angle)]]).dot(np.asarray((1.0, 0.0)))
        N = np.asarray(N)[0]
        solid.__init__(self, P, N, np.asarray((0.0, 0.0)), 0.0)
        
        self.__m = m
        self.__I = I
        self.__z = z_cog
        self.__p = p
        self.__h = h
        self.__w = w
        self.alpha = 0.0

    @property
    def p(self):
        """Propeller"""
        return self.__p

    @property
    def h(self):
        """Hull"""
        return self.__h

    @property
    def w(self):
        """Wings"""
        return self.__w

    def update(self, dt=DT):
        """Time integration.
        """
        # Get the current mass and center of gravity (they could be modified by the water)
        water = self.__p.water
        m_w = water.getM()
        z_w = water.getC()
        m = self.__m + m_w
        z = (self.__z * self.__m + z_w * m_w) / m
        
        # Modify the inertia due to the water mass and cog change
        z_0 = np.sqrt(self.__I / self.__m)
        dz = self.__z - z
        I = self.__m * (z_0 + dz)**2 + m_w * (z_w - z)**2

        # Get the attack angle
        U = np.linalg.norm(self.U)
        if U == 0.0:
            alpha = 0.0
        else:
            n = self.U.dot(self.N)
            t = self.U.dot(np.asarray((-self.N[1], self.N[0])))
            alpha = np.arctan2(t, n)
        self.alpha = alpha
        
        # Compute the corrected gravity
        gn = G * self.N[1]
        gt = -G * self.N[0]  # T = (N[1], -N[0])
        
        # Get the forces and its application points
        F = []
        z_F = []
        F.append(np.asarray((gt, gn)))
        z_F.append(z)
        F.append(np.asarray((0.0, __NU_PROP__ * self.__p.getFprop())))
        z_F.append(z)
        F.append(np.asarray((self.__h.getL(alpha, U),
                             self.__h.getD(alpha, U))))
        z_F.append(self.__h.getZ(alpha, U))
        for w in self.__w:
            F.append(np.asarray((w.getL(alpha, U),
                                 w.getD(alpha, U))))
            z_F.append(w.getZ(alpha, U))
            
        # Integrate
        dUdt = solid.update(self, m, I, z, F, z_F, dt)
        self.__p.update(dUdt.dot(self.N), gn, dt)


    

In [70]:

    

# Setup the propeller
z_FS = __Z_FS__
p = __P__
bottle_test = propeller(geom_z, geom_r, z_FS, p)

# Setup the hull
alpha = np.asarray((np.radians(0.0), np.radians(5.0), np.radians(10.0), np.radians(15.0)))
CL = np.asarray((0.0, 0.016, 0.029, 0.035))
CD = np.asarray((0.023, 0.028, 0.033, 0.041))
Cz = np.asarray((0.68, 0.91, 0.98, 1.1))
hull_test = hull(hull_z, hull_r, alpha, CL, CD, Cz)

z = (0.0, 0.01)
r = (__L_WING__, __L_WING__)
wing_z = np.asarray(z)
wing_r = np.asarray(r)
wings_test = []

# Setup the horizontal wings
alphao = np.asarray((np.radians(0.0), np.radians(5.0), np.radians(10.0), np.radians(15.0)))
CLo = np.asarray((0.0, 0.5438, 0.9067, 0.4365))
CDo = np.asarray((0.0101, 0.0141, 0.0270, 0.1754))

wings_test.append(wing(wing_z, wing_r, alphao, CLo, CDo))
wings_test.append(wing(wing_z, wing_r, alphao, CLo, CDo))

# And the vertical ones
wings_test.append(wing(wing_z, wing_r, alphao, 0.0 * CLo, CDo))
wings_test.append(wing(wing_z, wing_r, alphao, 0.0 * CLo, CDo))

# Generate the rocket
mass = 5E-2
L_cyl = 0.175
cog = 0.8 * L_cyl
I = mass * (L_cyl - cog)**2
rocket_test = rocket(__ANGLE__, mass, I, cog, bottle_test, hull_test, wings_test)

t = [0.0]
x = [rocket_test.P[0]]
z = [rocket_test.P[1]]
a = [0.0]
while rocket_test.P[1] > 0.0:
    if rocket_test.p.water.z > 0.0:
        dt = DT
    else:
        dt = DT_BIG
    rocket_test.update(DT)
    t.append(t[-1] + DT)
    x.append(rocket_test.P[0])
    z.append(rocket_test.P[1])
    a.append(rocket_test.alpha)
    if rocket_test.alpha > np.radians(10.0) or rocket_test.alpha < -np.radians(10.0):
        break
plot(x, z, 'k-')
plot(x, np.degrees(a), 'r-')
grid()


    

In [64]:

    

# We must perform a initial shot to can start (with a minimum mass)
mass = 1E-2
L_cyl = 0.175
cog = 0.8 * L_cyl
I = mass * (L_cyl - cog)**2

# Setup the propeller
z_FS = __Z_FS__
p = __P__
bottle_test = propeller(geom_z, geom_r, z_FS, p)

# Setup the hull
hull_test = hull(hull_z, hull_r, alpha, CL, CD, Cz)

# Setup the horizontal wings
wings_test = []
wings_test.append(wing(wing_z, wing_r, alpha, CLo, CDo))
wings_test.append(wing(wing_z, wing_r, alpha, CLo, CDo))

# And the vertical ones
wings_test.append(wing(wing_z, wing_r, alpha, 0.0 * CLo, CDo))
wings_test.append(wing(wing_z, wing_r, alpha, 0.0 * CLo, CDo))

# Generate the rocket
rocket_test = rocket(0.5 * np.pi, mass, I, cog, bottle_test, hull_test, wings_test)

t = [0.0]
z = [rocket_test.P[1]]
while rocket_test.P[1] > 0.0 and rocket_test.U[1] >= 0.0:
    if rocket_test.p.water.z > 0.0:
        dt = DT
    else:
        dt = DT_BIG
    rocket_test.update(DT)
    t.append(t[-1] + DT)
    z.append(rocket_test.P[1])
plot(t, z, 'k-')
grid()

# Store the results
mass_max = mass
z_max = z[-1]
mm = [mass_max]
zz = [z_max]


    

In [65]:

    

# Now we can start checking configurations with other masses
dmass = 1E-2

while True:
    clear_output()
    print('mass = {:.2g} kg  ->  Z = {:.3g} m'.format(mass_max, z_max))
    mass = mass_max + dmass
    L_cyl = 0.175
    cog = 0.8 * L_cyl
    I = mass * (L_cyl - cog)**2
    print('checking = {:.2g} kg...'.format(mass))
    sys.stdout.flush()

    # Setup the propeller
    bottle_test = propeller(geom_z, geom_r, z_FS, p)

    # Setup the hull
    hull_test = hull(hull_z, hull_r, alpha, CL, CD, Cz)

    # Setup the horizontal wings
    wings_test = []
    wings_test.append(wing(wing_z, wing_r, alpha, CLo, CDo))
    wings_test.append(wing(wing_z, wing_r, alpha, CLo, CDo))

    # And the vertical ones
    wings_test.append(wing(wing_z, wing_r, alpha, 0.0 * CLo, CDo))
    wings_test.append(wing(wing_z, wing_r, alpha, 0.0 * CLo, CDo))

    # Generate the rocket
    rocket_test = rocket(0.5 * np.pi, mass, I, cog, bottle_test, hull_test, wings_test)

    t = [0.0]
    z = [rocket_test.P[1]]
    while rocket_test.P[1] > 0.0 and rocket_test.U[1] >= 0.0:
        if rocket_test.p.water.z > 0.0:
            dt = DT
        else:
            dt = DT_BIG
        rocket_test.update(DT)
        t.append(t[-1] + DT)
        z.append(rocket_test.P[1])

    mm.append(mass)
    zz.append(z[-1])
        
    if z[-1] < z_max:
        print('FAIL (Z = {:.2g} m)'.format(z_max))
        break
    # Store the results
    mass_max = mass
    z_max = z[-1]

plot(mm, zz, 'k-')
xlabel(r'$m$ [kg]')
ylabel(r'$z$ [m]')
grid()


    

    




mass = 0.25 kg  ->  Z = 56.1 m
checking = 0.26 kg...
FAIL (Z = 56 m)






    

In [66]:

    

r_max = max(geom_r)
z_id = numpy.where(geom_r == r_max)[0][0]
max_w_wing = geom_z[z_id]


    

In [67]:

    

# Again we are performing an initial shot to can start (with a minimum wing)
w_wing = 0.005

# Setup the propeller
bottle_test = propeller(geom_z, geom_r, z_FS, p)

# Setup the hull
hull_test = hull(hull_z, hull_r, alpha, CL, CD, Cz)

z = (0.0, w_wing)
r = (__L_WING__, __L_WING__)
wing_z = np.asarray(z)
wing_r = np.asarray(r)
wings_test = []

# Setup the horizontal wings
wings_test.append(wing(wing_z, wing_r, alphao, CLo, CDo))
wings_test.append(wing(wing_z, wing_r, alphao, CLo, CDo))

# And the vertical ones
wings_test.append(wing(wing_z, wing_r, alphao, 0.0 * CLo, CDo))
wings_test.append(wing(wing_z, wing_r, alphao, 0.0 * CLo, CDo))

# Generate the rocket
rocket_test = rocket(__ANGLE__, mass, I, cog, bottle_test, hull_test, wings_test)

t = [0.0]
x = [rocket_test.P[0]]
z = [rocket_test.P[1]]
a = [0.0]
while rocket_test.P[1] > 0.0:
    if rocket_test.p.water.z > 0.0:
        dt = DT
    else:
        dt = DT_BIG
    rocket_test.update(DT)
    t.append(t[-1] + DT)
    x.append(rocket_test.P[0])
    z.append(rocket_test.P[1])
    a.append(rocket_test.alpha)
    if rocket_test.alpha > np.radians(10.0) or rocket_test.alpha < -np.radians(10.0):
        break
plot(x, z, 'k-')
plot(x, np.degrees(a), 'r-')
grid()

# Store the result
wing_max = w_wing
x_max = x[-1]
ww = [wing_max]
xx = [x_max]


    

In [68]:

    

# Now we can start checking configurations with other masses
dwing = 5E-3

while True:
    clear_output()
    print('wing = {:.3g} m  ->  X = {:.3g} m'.format(wing_max, x_max))
    w_wing = wing_max + dwing
    print('checking = {:.3g} m...'.format(w_wing))
    sys.stdout.flush()

    if w_wing > max_w_wing:
        print('FAIL (The wing cannot be bigger than {} m)'.format(max_w_wing))
        break
    
    # Setup the propeller
    bottle_test = propeller(geom_z, geom_r, z_FS, p)

    # Setup the hull
    hull_test = hull(hull_z, hull_r, alpha, CL, CD, Cz)

    z = (0.0, w_wing)
    r = (__L_WING__, __L_WING__)
    wing_z = np.asarray(z)
    wing_r = np.asarray(r)
    wings_test = []

    # Setup the horizontal wings
    wings_test.append(wing(wing_z, wing_r, alpha, CLo, CDo))
    wings_test.append(wing(wing_z, wing_r, alpha, CLo, CDo))

    # And the vertical ones
    wings_test.append(wing(wing_z, wing_r, alpha, 0.0 * CLo, CDo))
    wings_test.append(wing(wing_z, wing_r, alpha, 0.0 * CLo, CDo))

    # Generate the rocket
    rocket_test = rocket(__ANGLE__, mass, I, cog, bottle_test, hull_test, wings_test)

    t = [0.0]
    x = [rocket_test.P[0]]
    z = [rocket_test.P[1]]
    a = [0.0]
    alpha_valid = True
    while rocket_test.P[1] > 0.0:
        if rocket_test.p.water.z > 0.0:
            dt = DT
        else:
            dt = DT_BIG
        rocket_test.update(DT)
        t.append(t[-1] + DT)
        x.append(rocket_test.P[0])
        z.append(rocket_test.P[1])
        a.append(rocket_test.alpha)
        if rocket_test.alpha > np.radians(10.0) or rocket_test.alpha < -np.radians(10.0):
            alpha_valid = False
            break

    ww.append(w_wing)
    xx.append(x[-1])

    if x[-1] < x_max and alpha_valid:
        print('FAIL (X = {:.3g} m)'.format(x[-1]))
        break
    # Store the results
    wing_max = w_wing
    x_max = x[-1]

plot(ww, xx, 'k-')
xlabel(r'$wing$ [m]')
ylabel(r'$x$ [m]')
grid()


    

    




wing = 0.06 m  ->  X = 40.2 m
checking = 0.065 m...
FAIL (The wing cannot be bigger than 0.06 m)






    

In [69]:

    

plot(x, z, 'k-')
plot(x, np.degrees(a), 'r-')
xlabel(r'$x$ [m]')
ylabel(r'$z$ [m]')
grid()


    

In [ ]: