diff --git a/src/debug_mecanica_disparo.py b/src/debug_mecanica_disparo.py new file mode 100644 index 0000000..6f67cdf --- /dev/null +++ b/src/debug_mecanica_disparo.py @@ -0,0 +1,74 @@ +from funcionesSimELectrica import * +from funcionesSimFisica import * +from scipy.integrate import cumulative_trapezoid + +import matplotlib.pyplot as plt +import numpy as np + + +__author__ = 'Oscar Suescun' + +lts_path = "C:/Users/osuescuneli/AppData/Local/Programs/ADI/LTspice/LTspice.exe" + +''' +Esto sirve para hacer debug del sistema y ver como se comportan +las funciones para la + +''' +if __name__ == '__main__': + tiempoSimulacion = 10e-3 + step = tiempoSimulacion / 1e6 + + numeroModulos = 2 + + longitudBobina = 53.21e-3 # m + 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-3 # F + + masa = calcular_masa(diametroVastago, longitudVastago) + anguloDisparo = 40 + + + resistencia = calcular_resistencia(diametroCuBobina, longitudBobina, + diametroInteriorBobina, espirasBobina) + + reluctancia = np.mean(reluctancia_funcX(longitudVastago, diametroVastago, + longitudBobina, diametroInteriorBobina, + diametroCuBobina, espirasBobina)) + + inductancia = espirasBobina**2 / reluctancia + + tiempo = np.arange(0,tiempoSimulacion, step) + + corriente = solver_manual(Tension, resistencia, inductancia, Capacitancia, tiempo) + fuerza = fuerza_magnetica(corriente, reluctancia, espirasBobina, diametroVastago) + + aceleracion = [] + + for i in np.arange(numeroModulos): + + aceleracion = np.append(aceleracion, calcular_aceleracion(masa, anguloDisparo, fuerza)) + + if i == 0:aceleracion[aceleracion<0] = 0 # Evitar que la aceleracion sea negativa en la primera vuelta + else:pass + + velocidad = cumulative_trapezoid(aceleracion,tiempo,initial=0) + posicion = cumulative_trapezoid(velocidad, tiempo, initial=0) + + idx = posicion <= longitudBobina*(i+1) + + tiempo = tiempo[idx] # Recorto los vectores + posicion = posicion[idx] # para tener en cuenta + velocidad = velocidad[idx] # que el vastago recorre + aceleracion = aceleracion[idx] # la longitud de la bobina + + + + diff --git a/src/funcionesSim.py b/src/funcionesSim.py deleted file mode 100644 index 6f3deee..0000000 --- a/src/funcionesSim.py +++ /dev/null @@ -1,114 +0,0 @@ -from PyLTSpice import SimRunner, SpiceEditor, RawRead -import numpy as np -import matplotlib.pyplot as plt - -''' -Este codigo funciona sobre modelo.net - -Para conseguir este archivo hay que ir al .asc con -el mismo nombre y ir View > SPICE Netlist. Esto -mostrara el codigo que se quiere. Se copia el codigo -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): - - ''' - Funcion principal de simulacion. - - Entradas Obligatorias: - - tensionCap -> Tension del capacitor - - Rtotal -> Resistencia total en el circuito (Rextra + Rbobina) - - Lbobina -> Inductancia de la bobina - - Ctotal -> capacidad del banco de capacitores - - lts_path -> directorio donde tienes el .exe de LTSpice - - Entradas Opcionales: - - tSim -> tiempo total de simulacion (de normal 100ms) - - toff_sw -> tiempo que tarda la fuente en - desconectarse del capacitor (de normal 1ms) - - ton_mos -> tiempo que tarda el mosfet en abrir el gate - del mosfet de disparo (de normal 1.5ms) - - Salida - resultados = { - 'tiempo' : vector con el tiempo de simulacion - 'vCap' : Vector con la tension entre Capacitor -> GND - 'vBob' : Vector de tension atraves de la bobina - 'iBob' : Vector con la corriente atraves de Bobina - } - - ''' - - modeloSim = "simulador/modelo_transitorio.asc" - outDir = 'simulador' - - ## Hago la simulacion - - simulador = SimRunner(output_folder=outDir, simulator=lts_path) - - simulador.create_netlist(modeloSim) - netlist = SpiceEditor(modeloSim.replace('.asc', '.net')) - - netlist.set_parameters(V = tensionCap, - R = Rtotal, - L = Lbobina, - C = Ctotal, - tSim = tSim, - toff_sw = toff_sw, - ton_mos = ton_mos) - - simulador.run(netlist) - - raw_path, log_path = next(iter(simulador)) # Obtener el único resultado - raw = RawRead(raw_path) - - 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)')) - } - except KeyError as e: - print(f'Error al leer señal del archivo .raw: {e}') - return None - -def dibujar(resultado): - plt.figure() - plt.plot(resultado['tiempo'], resultado['vCap'], label='Tension Cap [V]') - plt.plot(resultado['tiempo'], resultado['vBob'], label='Tension Bobina [V]') - plt.plot(resultado['tiempo'], resultado['iBob'], label='Corriente Bobina [A]') - plt.grid() - plt.legend() - plt.show() - - -if __name__ == '__main__': - - lts_path = "D:\\Appdata\\LTSpice\\LTSPice.exe" - - Tension = 50 - Resistencia = 1.942 - Inductancia = 11e-6 - Capacitancia = 1e-3 - - tiempoSimulacion = 15e-3 - tiempoSwitch = 25e-6 - tiempoMos = 50e-6 - - - resultado = simular(Tension, - Resistencia, - Inductancia, - Capacitancia, - lts_path, - tSim = tiempoSimulacion, - toff_sw = tiempoSwitch, - ton_mos = tiempoMos) - - dibujar(resultado) - \ No newline at end of file diff --git a/src/funcionesSimELectrica.py b/src/funcionesSimELectrica.py new file mode 100644 index 0000000..9fbe9fb --- /dev/null +++ b/src/funcionesSimELectrica.py @@ -0,0 +1,169 @@ +from PyLTSpice import SimRunner, SpiceEditor, RawRead +import numpy as np +import matplotlib.pyplot as plt + +''' +Este codigo funciona sobre modelo.net + +Para conseguir este archivo hay que ir al .asc con +el mismo nombre y ir View > SPICE Netlist. Esto +mostrara el codigo que se quiere. Se copia el codigo +y se pega en un file con extenion .net + +''' + +__author__ = 'Oscar Suescun' + +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. Es legado por si se quisiera simular + con LTSpice. No se usa porque hay problemas con el tamaño del step + en los resultados. + + Entradas Obligatorias: + - tensionCap -> Tension del capacitor + - Rtotal -> Resistencia total en el circuito (Rextra + Rbobina) + - Lbobina -> Inductancia de la bobina + - Ctotal -> capacidad del banco de capacitores + - lts_path -> directorio donde tienes el .exe de LTSpice + + Entradas Opcionales: + - tSim -> tiempo total de simulacion (de normal 100ms) + - toff_sw -> tiempo que tarda la fuente en + desconectarse del capacitor (de normal 1ms) + - ton_mos -> tiempo que tarda el mosfet en abrir el gate + del mosfet de disparo (de normal 1.5ms) + + Salida + resultados = { + 'tiempo' : vector con el tiempo de simulacion + 'vCap' : Vector con la tension entre Capacitor -> GND + 'vBob' : Vector de tension atraves de la bobina + 'iBob' : Vector con la corriente atraves de Bobina + } + + ''' + + modeloSim = "simulador/modelo_transitorio.asc" + outDir = 'simulador' + + ## Hago la simulacion + + simulador = SimRunner(output_folder=outDir, simulator=lts_path) + + simulador.create_netlist(modeloSim) + netlist = SpiceEditor(modeloSim.replace('.asc', '.net')) + + netlist.set_parameters(V = tensionCap, + R = Rtotal, + L = Lbobina, + C = Ctotal, + tSim = tSim, + toff_sw = toff_sw, + ton_mos = ton_mos) + + simulador.run(netlist) + + 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': 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}') + return None + +def dibujar(resultado): + plt.figure() + plt.plot(resultado['tiempo']*1e3, resultado['vCap'], label='Tension Cap [V]') + plt.plot(resultado['tiempo']*1e3, resultado['vBob'], label='Tension Bobina [V]') + plt.plot(resultado['tiempo']*1e3, resultado['iBob'], label='Corriente Bobina [A]') + plt.xlabel("Tiempo [ms]") + plt.grid() + plt.legend() + plt.show() + + +if __name__ == '__main__': + + lts_path = "D:\\Appdata\\LTSpice\\LTSPice.exe" + + Tension = 30 + Resistencia = 10 + Inductancia = 1 + Capacitancia = 0.02 + + tiempoSimulacion = 10 + tiempoSwitch = 25e-6 + tiempoMos = 50e-6 + + + resultado1 = simular_LTS(Tension, + Resistencia, + Inductancia, + Capacitancia, + lts_path, + tSim = tiempoSimulacion, + toff_sw = tiempoSwitch, + ton_mos = tiempoMos) + + 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() + \ No newline at end of file diff --git a/src/funcionesSimFisica.py b/src/funcionesSimFisica.py new file mode 100644 index 0000000..879cde3 --- /dev/null +++ b/src/funcionesSimFisica.py @@ -0,0 +1,178 @@ +import numpy as np +import matplotlib.pyplot as plt + +global muVacio +global muFeRel +global kDisp +global cteRoz +global g +global resCu + +densidadHierro = 7800 # kg/m3 +muVacio = 4*np.pi*1e-7 +muFeRel = 5000 +resCu = 1.7e-8 # ohm/m +kDisp = 2 # cte de dispersion +cteRoz = 0.1 +g = 9.8 + +__author__ = 'Oscar Suescun' + +def calcular_diametro_ext(longitudBobina, diametroCuBobina, espirasBobina, diametroInteriorBobina): + + espirasCapa = int(longitudBobina/diametroCuBobina) + numeroCapas = int(espirasBobina/espirasCapa + .5) + + return 2 * numeroCapas * diametroCuBobina + diametroInteriorBobina + +def reluctancia_funcX(longFe, diamFe, longC, diamCint, diamCu, espiras): + ''' + Esta funcion devuelve un vector que contenga la reluctancia total en + funcion de la distancia que recorre el vastago dentro de la bobina. + + Para esto se utiliza un vector de 0 -> longC con un estep de longC/100 + + Entrada: + - longFe --> longitud del vastago + - diamFe --> diametro del vastago + - longC --> longitud de la bobina + - diamCint --> diametro interior de la bobina + - diamCu --> diametro del hilo de cobre + - espiras --> numero de espiras en la bobina + + Salida: + - relTotal --> devuelve un vector con la reluctancia total en + funcion de la distancia recorrida en la bobina. + + ''' + + diamCext = calcular_diametro_ext(longC, diamCu, espiras, diamCint) + + seccionFe = (np.pi * diamFe**2) / 4 + seccionSc = (np.pi * diamCext**2) / 4 + seccionDisp = (np.pi * (kDisp * diamCint)**2) / 4 + + x = np.arange(0,longC + longC/100, longC/100) + longFe = np.full(len(x), longFe) + longC = np.full(len(x), longC) + + + relFe = longFe / (muVacio * muFeRel * seccionFe) + relDisp = longC / (muVacio * seccionDisp) + relAire = (longC - x) / (muVacio * seccionSc) + relFlujo = (longC + longFe - x) / (muVacio * seccionDisp) + + relTotal = relFe + relDisp + relAire + relFlujo + + return relTotal + + +def induccion(corriente, reluctancia, espiras, diamFe): + flujo = espiras * corriente / reluctancia + + secFe = np.pi * diamFe**2 / 4 + + return flujo / secFe + +def fuerza_magnetica(corriente, reluctancia, espiras, diamFe): + ''' + Solo se calcula la fuerza magetica + + ''' + + b = induccion(corriente, reluctancia, espiras, diamFe) + + secFe = np.pi * diamFe**2 / 4 + + return np.sign(corriente) * 0.5 * b**2 * secFe / muVacio + +def calcular_aceleracion(masa, anguloDisparo, fuerzaMagnetica): + ''' + Retocar esto, tiene en cuenta el principio solo pero hay cosas raras + + ''' + + anguloDisparo = anguloDisparo * np.pi / 180 + acel = fuerzaMagnetica/masa - g*(cteRoz * np.cos(anguloDisparo) + np.sin(anguloDisparo)) + return acel + + +def calcular_resistencia(diamCu, longC, diamC, espiras): + + espirasCapa = int(longC/diamCu) + nCapasEnteras = int(espiras/espirasCapa) # solo capas completametne llenas + espirasUltimaCapa = espiras - espirasCapa * nCapasEnteras # las espiras de la ultima capa + + longitud = 0 + + for i in range(nCapasEnteras): longitud += espirasCapa*np.pi*(diamC + (1+i)*diamCu) # Calculo longitud de total de las capas llenas + + longitud += espirasUltimaCapa*np.pi*(diamC + (1 + nCapasEnteras)*diamCu) + + seccionCu = np.pi * diamCu**2 / 4 + + return resCu * longitud / seccionCu + +def calcular_masa(diamFe, longFe): return np.pi * diamFe**2 / 4 * longFe * densidadHierro + + + +if __name__ == '__main__': + from funcionesSimELectrica import * + from scipy.integrate import cumulative_trapezoid + + tiempoSimulacion = 10e-3 + step = tiempoSimulacion / 1e6 + + numeroModulos = 2 + + longitudBobina = 53.21e-3 # m + 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-3 # F + + masa = calcular_masa(diametroVastago, longitudVastago) + anguloDisparo = 40 + + + resistencia = calcular_resistencia(diametroCuBobina, longitudBobina, + diametroInteriorBobina, espirasBobina) + + reluctancia = np.mean(reluctancia_funcX(longitudVastago, diametroVastago, + longitudBobina, diametroInteriorBobina, + diametroCuBobina, espirasBobina)) + + inductancia = espirasBobina**2 / reluctancia + + tiempo = np.arange(0,tiempoSimulacion, step) + + corriente = solver_manual(Tension, resistencia, inductancia, Capacitancia, tiempo) + fuerza = fuerza_magnetica(corriente, reluctancia, espirasBobina, diametroVastago) + aceleracion = calcular_aceleracion(masa, anguloDisparo, fuerza) + aceleracion[aceleracion<0]=0 + velocidad = cumulative_trapezoid(aceleracion, tiempo, initial=0) + posicion = cumulative_trapezoid(velocidad, tiempo, initial=0) + + idx = posicion <= longitudBobina + + tiempo = tiempo[idx] + posicion = posicion[idx] + velocidad = velocidad[idx] + aceleracion = aceleracion[idx] + + plt.figure() + #plt.plot(t, corriente, label='corriente') + #plt.plot(t, fuerza, label = 'Fuerza') + plt.plot(tiempo, aceleracion, label='aceleracion') + # plt.plot(tiempo, velocidad, label = 'velocidad') + #plt.plot(posicion, velocidad, label = 'posicion') + plt.legend() + plt.show() + + diff --git a/src/simulador/modelo_transitorio.asc b/src/simulador/modelo_transitorio.asc index 2e2dfa1..453ed14 100644 --- a/src/simulador/modelo_transitorio.asc +++ b/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 diff --git a/src/simulador/modelo_transitorio.op.raw b/src/simulador/modelo_transitorio.op.raw new file mode 100644 index 0000000..9202c3f Binary files /dev/null and b/src/simulador/modelo_transitorio.op.raw differ diff --git a/src/simulador/modelo_transitorio.raw b/src/simulador/modelo_transitorio.raw new file mode 100644 index 0000000..2f8e36c Binary files /dev/null and b/src/simulador/modelo_transitorio.raw differ diff --git a/src/simulador/modelo_transitorio_1.net b/src/simulador/modelo_transitorio_1.net index ddaace6..eca5ccb 100644 --- a/src/simulador/modelo_transitorio_1.net +++ b/src/simulador/modelo_transitorio_1.net @@ -1,5 +1,10 @@ +<<<<<<< HEAD * C:\Users\pedro\Desktop\Projects\LaunchSim\src\simulador\modelo_transitorio.asc * Generated by LTspice 24.1.4 for Windows. +======= +* C:\Users\osuescuneli\Desktop\practicas\Practia_Lanzadera\source\src\simulador\modelo_transitorio.asc +* Generated by LTspice 24.1.5 for Windows. +>>>>>>> 04076e92921df836051a54c4429533fbd2cea3a7 V1 N001 0 {V} S1 N001 condensador cont1 0 SW V2 cont1 0 PULSE(0 10 0 1n 1n {toff_sw} {tsim}) @@ -8,18 +13,20 @@ 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\pedro\AppData\Local\LTspice\lib\cmp\standard.mos .tran 0 {tsim} 0 .Model SW SW(Ron=1m Roff=100Meg Vt = 5) -.param tsim 15m +.param tsim 10 .param toff_sw 25u .param ton_mos 50u -.param R 1.942 -.param C 1m -.param L 11u -.param V 50 +.param R 10 +.param C 20m +.param L 1 +.param V 30 .param delta tsim-ton_mos .backanno .end diff --git a/test.py b/test.py new file mode 100644 index 0000000..89feec2 --- /dev/null +++ b/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)