solver analitico añadido

src/debug_mecanica_disparo.py

@@ -16,22 +16,21 @@ las funciones para la
 
 '''
 if __name__ == '__main__':
-    tiempoSimulacion = 20e-3
-    tiempoSwitch = 1e-6
-    tiempoMos = 1e-3
+    tiempoSimulacion = 10e-3
+    step = tiempoSimulacion / 1e6
 
-    numeroModulos = 1
+    numeroModulos = 2
 
     longitudBobina = 53.21e-3 # m
-    diametroInteriorBobina = 6.035-3 * 2 # m
-    espirasBobina = 800
+    diametroInteriorBobina = 6.035e-3 * 2 # m
+    espirasBobina = 500
     diametroCuBobina = 0.8e-3 # m
 
     longitudVastago = 96e-3 # m
     diametroVastago = 3.045e-3 * 2 # m
         
     Tension = 30 # V
-    Capacitancia = 500e-6 # F
+    Capacitancia = 500e-3 # F
 
     masa = calcular_masa(diametroVastago, longitudVastago)
     anguloDisparo = 40
@@ -41,33 +40,46 @@ if __name__ == '__main__':
                                        diametroInteriorBobina, espirasBobina)
     
     reluctancia = np.mean(reluctancia_funcX(longitudVastago, diametroVastago,
-                                            longitudBobina, diametroInteriorBobina,
-                                            diametroCuBobina, espirasBobina))
+                                    longitudBobina, diametroInteriorBobina,
+                                    diametroCuBobina, espirasBobina))
     
     inductancia = espirasBobina**2 / reluctancia
 
-    resultadosSim = simular(Tension, resistencia,
-                            inductancia, Capacitancia, lts_path, 
-                            tSim= tiempoSimulacion,
-                            toff_sw= tiempoSwitch,
-                            ton_mos= tiempoMos)
-    
-    tiempo = resultadosSim['tiempo'][resultadosSim['tiempo'] >= tiempoMos]
+    tTemp = 0
+    tempAcel = 0
+    tempVelocidad = 0
+    longTemp = longitudBobina
 
-    dt = np.full(len(tiempo), tiempoMos)
-    tiempo -= dt
+    for i in np.arange(numeroModulos):
+        t = np.arange(0,tiempoSimulacion, step)
 
-    corriente = resultadosSim['iBob'][resultadosSim['tiempo'] >= tiempoMos]
+        corriente = solver_manual(Tension, resistencia, inductancia, Capacitancia, t)
+        fuerza = fuerza_magnetica(corriente, reluctancia, espirasBobina, diametroVastago)
+        acel = aceleracion(masa, 40, fuerza)
+
+        t = tTemp + t
+
+        velocidad = integrar(acel, t, initial = 0)
+        posicion = integrar(velocidad, t, initial = 0)
+
+        idx = posicion<=longTemp
+
+        t = t[idx]
+        corriente = corriente[idx]
+        fuerza = fuerza[idx]
+        acel = acel[idx]
+        velocidad = velocidad[idx]
+        posicion = posicion[idx]
+
+        tempAcel = acel[-1]
+        tempVelocidad = velocidad[-1]
+        longTemp = longitudBobina + longTemp
 
-    fuerzaMag = fuerza_magnetica(corriente, reluctancia, 
-                                 espirasBobina, diametroVastago)
-    
-    aVastago = aceleracion(masa, anguloDisparo, fuerzaMag)
-    vVastago = integrar(aVastago, initial=0)
-    xVastago = integrar(vVastago, initial=0)
 
     plt.figure()
-    plt.plot(tiempo * 1e3, xVastago)
-    plt.grid()
+    #plt.plot(t, corriente, label='corriente')
+    #plt.plot(t, fuerza, label = 'Fuerza')
+    #plt.plot(t, acel, label='aceleracion')
+    plt.plot(posicion, velocidad, label = 'velocidad')
     plt.legend()
     plt.show()

src/funcionesSimELectrica.py

@@ -14,7 +14,42 @@ y se pega en un file con extenion .net
 
 __author__ = 'Oscar Suescun'
 
-def simular(tensionCap, Rtotal, Lbobina, Ctotal, lts_path, tSim = 100e-3, toff_sw = 1e-3, ton_mos = 1.5e-3):
+def solver_manual(tensionCap, Rtotal, Lbobina, Ctotal, t):
+    '''
+    Solver analitico de la descarga LCR, resuelto manualmente.
+    Se esta resuelto con el metodo para ecuaciones diferenciales 
+    homogeneas de segundo orden con coeficientes constantes.
+
+    Para evitar que el sistema tenga un periodo oscilatorio se necesita
+    cumplir que R^2 >= 4 * L/C. Si esto no se cumple tendra un perido
+    de oscilacion antes de apagarse lo que puede generar fuerzas negativas
+    en el vastago.
+    
+    Entradas:
+        - tensionCap -> Tension del capacitor
+        - Rtotal -> Resistencia total en el circuito (Rextra + Rbobina)
+        - Lbobina -> Inductancia de la bobina
+        - Ctotal -> capacidad del banco de capacitores
+
+    Salida:
+        - Corriente en la bobina
+
+    '''
+
+    if Rtotal**2 - 4*Lbobina/Ctotal >= 0:
+        tau1 = (-Rtotal + np.sqrt(Rtotal**2 - 4*Lbobina/Ctotal)) / (2 * Lbobina)
+        tau2 = (-Rtotal - np.sqrt(Rtotal**2 - 4*Lbobina/Ctotal)) / (2 * Lbobina)
+        iCap = tensionCap/(Lbobina*(tau1 - tau2)) * (np.exp(tau1 * t) - np.exp(tau2 * t))
+
+    else:
+        a = -Rtotal / (2*Lbobina)
+        b = np.sqrt(4*Lbobina/Ctotal - Rtotal**2) / (2*Lbobina)
+        iCap = np.exp(a*t) * tensionCap / (b * Lbobina) * np.sin(b*t)
+
+    return iCap
+
+
+def simular_LTS(tensionCap, Rtotal, Lbobina, Ctotal, lts_path, tSim = 100e-3, toff_sw = 1e-3, ton_mos = 1.5e-3):
 
     '''
     Funcion principal de simulacion.
@@ -66,13 +101,21 @@ def simular(tensionCap, Rtotal, Lbobina, Ctotal, lts_path, tSim = 100e-3, toff_s
     raw_path, log_path = next(iter(simulador))  # Obtener el único resultado
     raw = RawRead(raw_path)
 
+    tiempo = np.array(raw.get_trace('time'))
+    vCap = np.array(raw.get_trace(('V(condensador)')))
+    vBob = np.array(raw.get_trace('V(v1)')) - np.array(raw.get_trace('V(v2)'))
+    iBob = np.array(raw.get_trace('I(L1)'))
+    iMos = np.array(raw.get_trace('Id(M1)'))
+
+    idx = np.where(tiempo >= ton_mos)
+
     try:
         return {
-            'tiempo': np.array(raw.get_trace('time')),
-            'vCap': np.array(raw.get_trace('V(condensador)')),
-            'vBob': np.array(raw.get_trace('V(v1)')) - np.array(raw.get_trace('V(v2)')),
-            'iBob': np.array(raw.get_trace('I(L1)')),
-            'iMos' : np.array(raw.get_trace('Id(M1)'))
+            'tiempo': tiempo[idx],
+            'vCap': vCap[idx],
+            'vBob': vBob[idx],
+            'iBob': iBob[idx],
+            'iMos' : iMos[idx]
         }
     except KeyError as e:
         print(f'Error al leer señal del archivo .raw: {e}')
@@ -93,17 +136,17 @@ if __name__ == '__main__':
 
     lts_path = "C:/Users/osuescuneli/AppData/Local/Programs/ADI/LTspice/LTspice.exe"
 
-    Tension = 150
-    Resistencia = 4
-    Inductancia = 100e-6
-    Capacitancia = 500e-6
+    Tension = 30
+    Resistencia = 10
+    Inductancia = 1
+    Capacitancia = 0.02
 
-    tiempoSimulacion = 15e-3
+    tiempoSimulacion = 10
     tiempoSwitch = 25e-6
     tiempoMos = 50e-6
 
 
-    resultado = simular(Tension, 
+    resultado1 = simular_LTS(Tension, 
                         Resistencia, 
                         Inductancia, 
                         Capacitancia, 
@@ -112,5 +155,13 @@ if __name__ == '__main__':
                         toff_sw = tiempoSwitch,
                         ton_mos = tiempoMos)
     
-    dibujar(resultado)
+    t = np.arange(0,tiempoSimulacion, 10e-6)
+
+    resultado2 = solver_manual(Tension, Resistencia, Inductancia, Capacitancia, t )
+    
+    plt.figure()
+    plt.plot(resultado1['tiempo']*1e3, resultado1['iBob'], label = 'LTS')
+    plt.plot(t*1e3, resultado2, label='Manual')
+    plt.legend()
+    plt.show()
     

src/funcionesSimFisica.py

@@ -21,7 +21,7 @@ __author__ = 'Oscar Suescun'
 def calcular_diametro_ext(longitudBobina, diametroCuBobina, espirasBobina, diametroInteriorBobina):
     
     espirasCapa = int(longitudBobina/diametroCuBobina)
-    numeroCapas = int(espirasBobina/espirasCapa +.5)
+    numeroCapas = int(espirasBobina/espirasCapa + .5)
 
     return 2 * numeroCapas * diametroCuBobina + diametroInteriorBobina
 
@@ -75,6 +75,10 @@ def induccion(corriente, reluctancia, espiras, diamFe):
     return flujo / secFe
 
 def fuerza_magnetica(corriente, reluctancia, espiras, diamFe):
+    '''
+    Solo se calcula la fuerza magetica
+
+    '''
 
     b = induccion(corriente, reluctancia, espiras, diamFe)
 
@@ -83,10 +87,16 @@ def fuerza_magnetica(corriente, reluctancia, espiras, diamFe):
     return  np.sign(corriente) * 0.5 * b**2 * secFe / muVacio
 
 def aceleracion(masa, anguloDisparo, fuerzaMagnetica):
+    '''
+    Retocar esto, tiene en cuenta el principio solo pero hay cosas raras
+    
+    '''
 
     anguloDisparo = anguloDisparo * np.pi / 180
 
-    return fuerzaMagnetica/masa - g*(cteRoz * np.cos(anguloDisparo) + np.sin(anguloDisparo))
+    acel = fuerzaMagnetica/masa - g*(cteRoz * np.cos(anguloDisparo) + np.sin(anguloDisparo))
+    acel[acel < 0 ] = 0
+    return acel
     
     
 def calcular_resistencia(diamCu, longC, diamC, espiras):
@@ -113,7 +123,7 @@ if __name__ == '__main__':
     numeroModulos = 1
 
     longitudBobina = 53.21e-3 # m
-    diametroInteriorBobina = 6.035-3 * 2 # m
+    diametroInteriorBobina = 6.035e-3 * 2 # m
     espirasBobina = 600
     diametroCuBobina = 0.8e-3 # m
 

src/simulador/modelo_transitorio.asc

@@ -13,6 +13,12 @@ WIRE -256 240 -256 160
 WIRE 256 304 256 272
 WIRE 256 320 256 304
 WIRE -256 336 -256 320
+WIRE 208 608 64 608
+WIRE 64 672 64 608
+WIRE 208 672 208 608
+WIRE 64 800 64 752
+WIRE 208 800 208 752
+WIRE 208 800 64 800
 FLAG 432 192 0
 FLAG 432 112 cont1
 FLAG -208 96 cont1
@@ -59,13 +65,21 @@ SYMATTR Value {L}
 SYMBOL res 240 64 R0
 SYMATTR InstName R1
 SYMATTR Value {R}
-TEXT -320 -16 Left 2 !.tran 0 {tsim} 0
+SYMBOL voltage 64 656 R0
+WINDOW 123 0 0 Left 0
+WINDOW 39 0 0 Left 0
+SYMATTR InstName V4
+SYMATTR Value 30
+SYMBOL res 192 656 R0
+SYMATTR InstName R2
+SYMATTR Value 20
+TEXT -320 -16 Left 2 !.tran 0 {tsim} 0 
 TEXT -320 -48 Left 2 !.Model SW SW(Ron=1m Roff=100Meg Vt = 5)
-TEXT -320 -176 Left 2 !.param tsim 100m
+TEXT -320 -176 Left 2 !.param tsim 1.6m
 TEXT -320 -144 Left 2 !.param toff_sw 1m
 TEXT -320 -112 Left 2 !.param ton_mos 1.5m
-TEXT 336 -96 Left 2 !.param R 1.942
-TEXT 336 -160 Left 2 !.param C 1m
-TEXT 336 -128 Left 2 !.param L 11u
-TEXT 336 -64 Left 2 !.param V 150
+TEXT 336 -96 Left 2 !.param R 10
+TEXT 336 -160 Left 2 !.param C 0.02
+TEXT 336 -128 Left 2 !.param L 1
+TEXT 336 -64 Left 2 !.param V 30
 TEXT -320 -80 Left 2 !.param delta tsim-ton_mos

src/simulador/modelo_transitorio.net

@@ -1,25 +0,0 @@
-* C:\Users\osuescuneli\Desktop\source\src\simulador\modelo_transitorio.asc
-* Generated by LTspice 24.1.5 for Windows.
-V1 N001 0 {V}
-S1 N001 condensador cont1 0 SW
-V2 cont1 0 PULSE(0 10 0 1n 1n {toff_sw} {tsim})
-C1 condensador 0 {C}
-M1 V2 cont2 0 0 IRFZ44N
-V3 cont2 0 PULSE(0 15 {ton_mos} 1n 1n {delta} {tsim})
-L1 V1 V2 {L}
-R1 condensador V1 {R}
-.model NMOS NMOS
-.model PMOS PMOS
-.lib C:\Users\osuescuneli\AppData\Local\LTspice\lib\cmp\standard.mos
-.tran 0 {tsim} 0
-.Model SW SW(Ron=1m Roff=100Meg Vt = 5)
-.param tsim 100m
-.param toff_sw 1m
-.param ton_mos 1.5m
-.param R 1.942
-.param C 1m
-.param L 11u
-.param V 150
-.param delta tsim-ton_mos
-.backanno
-.end

src/simulador/modelo_transitorio_1.net

@@ -1,4 +1,4 @@
-* C:\Users\osuescuneli\Desktop\source\src\simulador\modelo_transitorio.asc
+* C:\Users\osuescuneli\Desktop\practicas\Practia_Lanzadera\source\src\simulador\modelo_transitorio.asc
 * Generated by LTspice 24.1.5 for Windows.
 V1 N001 0 {V}
 S1 N001 condensador cont1 0 SW
@@ -8,17 +8,19 @@ M1 V2 cont2 0 0 IRFZ44N
 V3 cont2 0 PULSE(0 15 {ton_mos} 1n 1n {delta} {tsim})
 L1 V1 V2 {L}
 R1 condensador V1 {R}
+V4 N002 N003 30
+R2 N002 N003 20
 .model NMOS NMOS
 .model PMOS PMOS
 .lib C:\Users\osuescuneli\AppData\Local\LTspice\lib\cmp\standard.mos
 .tran 0 {tsim} 0
 .Model SW SW(Ron=1m Roff=100Meg Vt = 5)
-.param tsim 20m
-.param toff_sw 1u
-.param ton_mos 1m
-.param R 3.42142
-.param C 500u
-.param L 13.8963m
+.param tsim 10
+.param toff_sw 25u
+.param ton_mos 50u
+.param R 10
+.param C 20m
+.param L 1
 .param V 30
 .param delta tsim-ton_mos
 .backanno

test.py

@@ -0,0 +1,5 @@
+import numpy as np
+
+vector = np.array([-3, -1, 0, 2, 5, -7])
+vector[vector < 0] = 0
+print(vector)