Aprender › Ingeniería

Simulación multidominio: cómo calcula el simulador

El simulador multidominio junta en un mismo esquema circuitos, motores, calefactores, tanques, sensores y bloques de control, y los resuelve todos a la vez en el tiempo. No tiene resolución detallada en PDF porque no hay una «solución a mano»: es un integrador numérico. Este tema explica qué hace por dentro, para que sepas leer sus resultados, confiar en ellos y reconocer sus límites.

Ir al simulador

Un mismo lenguaje para varios dominios

Todos los sistemas físicos que simula se describen con dos tipos de variables: una de esfuerzo o diferencia (voltaje, velocidad angular, temperatura, nivel o presión) y una de flujo (corriente, par, flujo de calor, caudal). Con ellas, los elementos de cada dominio se parecen a los del circuito eléctrico, y sus ecuaciones tienen la misma forma:

Dominio	Esfuerzo · flujo	Resistencia	Almacena (capacidad)	Ecuación del almacén
Eléctrico	V · i	$R$ : $V = R i$	capacitor $C$	$i = C \frac{d V}{d t}$
Mecánico (giro)	ω · T	fricción $b$ : $T = bω$	inercia $J$	$T = J \frac{d ω}{d t}$
Térmico	T · q	$R_{t h}$ : $Δ T = R_{t h} q$	capacidad $C_{t h}$	$q = C_{t h} \frac{d T}{d t}$
Hidráulico	h · Q	salida $Q = k h$ (no lineal)	área del tanque $A$	$Q = A \frac{d h}{d t}$

Por eso un circuito RC, un calefactor que se calienta y un tanque que se llena responden igual, con una exponencial de constante de tiempo $τ = R \cdot C$ (o $R_{t h} C_{t h}$ , o $2 A h / k$ ). Los componentes que convierten energía unen dos dominios: el motor de corriente directa produce un par $T_{m} = K i$ y a la vez genera una fuerza contraelectromotriz $K ω$ :

L \frac{d i}{d t} + R i + K ω = v J \frac{d ω}{d t} + b ω + T_{L} = K i

Motor de corriente directa: la corriente da par, la velocidad da fuerza contraelectromotriz, y las dos ecuaciones se resuelven juntas.

Cómo arma el sistema de ecuaciones

Las incógnitas son los voltajes de los nodos del esquema y las variables de estado de los componentes que guardan energía (el voltaje de un capacitor, la corriente de una bobina, la corriente, la velocidad y el ángulo de un motor, la temperatura de un calefactor, el nivel de un tanque…). Todas van en un vector $z$ . Cada componente aporta dos cosas:

Las corrientes que entran o salen por sus terminales. En cada nodo la suma tiene que ser cero (ley de corrientes de Kirchhoff): es el análisis nodal.
Sus ecuaciones internas, escritas como «residuo = 0»: para un capacitor $i - C \overset{v}{˙} = 0$ ; para una bobina $L \overset{}{˙} + R_{s} i - v = 0$ ; para un calefactor $C_{t h} \dot{T} - [P - (T - T_{a}) / R_{t h}] = 0$ .

Juntando todo queda un sistema de ecuaciones algebraicas y diferenciales con la forma general

F (z, \overset{z}{˙}, t) = 0

Este formato admite cualquier componente sin tratarlo aparte: un diodo (ecuación de Shockley), un transistor, un amplificador operacional con saturación o un sensor con su propio retraso son solo más ecuaciones en $F$ . Las fuentes de voltaje tienen una resistencia interna pequeña (0.01 Ω por omisión), como las reales, y se añade una conductancia diminuta de cada nodo a tierra ( $1 0^{- 12}$ S) para que un nodo aislado no deje el sistema sin solución.

Integración en el tiempo: BDF2 y Newton

El tiempo avanza en pasos $h$ . En cada paso las derivadas $\overset{z}{˙}$ se sustituyen por una fórmula con los valores de los pasos anteriores, y lo que queda es un sistema algebraico en $z_{n + 1}$ . El simulador usa la fórmula de diferencias hacia atrás de orden 2 (BDF2), que toma la pendiente en el punto nuevo de la parábola que pasa por los tres últimos:

\overset{z}{˙}_{n + 1} \approx \frac{3 z _{n + 1} - 4 z _{n} + z _{n - 1}}{2 h}

BDF2: la derivada en tₙ₊₁ es la pendiente de la parábola por los tres últimos puntos.

BDF2 es implícito: la incógnita aparece en los dos lados, así que hay que resolver $F (z_{n + 1}) = 0$ en cada paso. Se hace con el método de Newton: se evalúa el residuo $F$ , se calcula su jacobiano $J = \partial F / \partial z$ (derivando numéricamente cada componente), se resuelve $J Δ z = - F$ y se corrige $z$ hasta que la corrección es despreciable. Los métodos implícitos son estables aunque el sistema tenga dinámicas muy rápidas (un circuito con constantes de microsegundos dentro de un proceso térmico de minutos): el paso no tiene que ser más pequeño que la constante más rápida.

Un paso del simulador. El primer paso (sin historia) usa Euler hacia atrás; si Newton no converge, el paso se divide en dos mitades, hasta 8 veces.

El tamaño del paso. Por omisión es la duración de la simulación entre 2000, y se reduce si hay fuentes rápidas, muestreos o retardos. Además se ajusta a las oscilaciones del circuito: BDF2 amortigua una oscilación de frecuencia $ω$ en un factor de aproximadamente $(ωh)^{4} /4$ por paso y la retrasa unos $(2/9) (ωh)^{3}$ radianes, así que el simulador elige $h$ para que, en toda la simulación, ese amortiguamiento numérico quede por debajo del 1 % y el desfase por debajo de 0.01 rad. Los componentes digitales (relojes, contadores, el 555, comparadores con histéresis) cambian por eventos al final de cada paso aceptado.

Linealización: polos, estabilidad y G(s)

Al final de la simulación (o en el instante elegido) el simulador linealiza el sistema alrededor del estado en que quedó: deriva numéricamente las ecuaciones respecto a los estados y a la entrada, y obtiene un modelo lineal

\overset{x}{˙} = A x + B u y = C x + D u G (s) = C (s I - A)^{- 1} B + D

Los valores propios de A son los polos del sistema completo; si todos tienen parte real negativa, ese punto de operación es estable. Con la misma matriz se calculan los polos propios de cada subsistema (la red de capacitores y bobinas, un motor, un PID…), la función de transferencia entre la entrada y la salida que elijas, su respuesta al escalón, el diagrama de Bode y el lugar de las raíces respecto a cualquier parámetro (se vuelve a linealizar para cada valor).

Un modelo lineal solo vale cerca del punto donde se linealizó. Un diodo, un motor saturado, un control ON/OFF o la salida

k h

de un tanque cambian de comportamiento lejos de ese punto: por eso el simulador compara, cuando puede, la respuesta del modelo lineal con la simulación no lineal, y avisa cuando los polos dicen «estable» pero la simulación oscila (un ciclo límite).

Cómo leer y comprobar los resultados

Unidades y orden de magnitud: antes de mirar decimales, comprueba con una cuenta rápida (V/R, $K V / (R b + K^{2})$ , $P R_{t h}$ ) que el resultado tiene sentido.
Duración: simula al menos 5 constantes de tiempo de la dinámica más lenta para ver el valor final; la tarjeta de polos te da esas constantes ( $τ = 1/∣ σ ∣$ ).
Paso: si dudas de un detalle fino (un pico, un rizado), repite con un paso menor; si el resultado no cambia, es del modelo y no del método numérico.
Avisos: el simulador avisa de componentes fuera de sus límites (una resistencia que se quemaría, un motor bloqueado, un tanque que se desborda) y de esquemas mal planteados (fuentes en cortocircuito o en paralelo con valores distintos).

Errores frecuentes

Olvidar la tierra: sin un nodo de referencia los voltajes no están definidos.
Leer un valor «final» antes de que el sistema se asiente: mira la respuesta completa o alarga la simulación.
Usar la función de transferencia lejos del punto de operación de un sistema no lineal.
Esperar la precisión de un cálculo exacto en el instante de una conmutación: justo después de un salto de la fuente, BDF2 mezcla la historia de antes del salto durante un par de pasos.
Confundir el polo de un integrador (en 0, «marginal») con una inestabilidad: lo que importa es el lazo cerrado.

Ejemplos resueltos

Ejemplo 1Carga y descarga de un capacitor

Una fuente de pulsos de 0 a 5 V y 2 Hz alimenta una resistencia de 10 kΩ en serie con un capacitor de 10 µF. Halla la constante de tiempo, el polo y los voltajes máximo y mínimo del capacitor en régimen periódico.

Solución

$τ = RC = 1 0^{4} \cdot 1 0^{- 5} = 0.1$ s; el polo está en $- 1/ τ = - 10$ s⁻¹ y $G (s) = V_{C} / V = \frac{10}{s + 10}$ .
Cada medio periodo dura 0.25 s = 2.5τ. En régimen periódico, la carga va de $V_{min}$ a $V_{ma x}$ y la descarga vuelve: $V_{ma x} = 5 + (V_{min} - 5) e^{- 2.5}$ y $V_{min} = V_{ma x} e^{- 2.5}$ .
Resolviendo: $V_{ma x} = \frac{5 ( 1 - e ^{- 2.5} )}{1 - e ^{- 5}} = 4.621$ V y $V_{min} = 4.621 e^{- 2.5} = 0.379$ V.
Respuesta del modelo lineal a un escalón: subida 10–90 % en $τ ln 9 = 220$ ms y establecimiento al 2 % en $τ ln 50 = 391$ ms.

τ = 0.1

s, polo

- 10

; en régimen periódico

V_{C}

oscila entre 0.379 y 4.621 V.

Resuélvelo en Simuveliq

En Ingeniería → Simulación multidominio, en Ejemplos, «Cargar ejemplo» → Carga y descarga de un capacitor (RC). El Osciloscopio muestra la onda de la fuente y la del capacitor (máximos de 4.62 V y mínimos de 0.38 V con el paso de 1 ms).
Estabilidad: sistema completo y subsistemas da el polo −10 con τ = 100 ms y «Sigue a la fuente periódica»; Análisis de control da G(s) = 10/(s + 10) y, para el modelo lineal, subida 220 ms, establecimiento 392 ms y valor final 1.

Abrir en el simulador →

Ejemplo 2Circuito RLC subamortiguado

Un pulso de 1 V alimenta R = 10 Ω, L = 10 mH y C = 100 µF en serie (la fuente tiene 0.01 Ω internos). Halla los polos, el amortiguamiento, el sobreimpulso del voltaje del capacitor y la frecuencia de la oscilación.

Solución

Malla: $L \overset{q}{¨} + R \overset{q}{˙} + q / C = v$ , así que $\frac{V _{C}}{V} = \frac{1/ L C}{s ^{2} + ( R / L ) s + 1/ L C}$ .
$ω_{n} = 1/ L C = 1000$ rad/s. Con $R = 10.01$ Ω (10 + 0.01 de la fuente): $σ = R /2 L = 500.5$ s⁻¹, $ζ = σ / ω_{n} = 0.5005$ y $ω_{d} = ω_{n}^{2} - σ^{2} = 865.7$ rad/s. Polos: $- 500.5 \pm 865.7 j$ .
$M_{p} = e^{- π ζ / 1 - ζ^{2}} = 16.3 %$ : el capacitor llega a 1.163 V; $t_{p} = π / ω_{d} = 3.63$ ms. La oscilación tiene $f_{d} = ω_{d} /2 π = 137.8$ Hz.

Polos

- 500.5 \pm 865.7 j

(

ζ = 0.5

ω_{n} = 1000

rad/s);

V_{C}

máximo 1.163 V a los 3.6 ms del flanco.

Resuélvelo en Simuveliq

«Cargar ejemplo» → Circuito RLC subamortiguado (polos complejos). El polo dominante sale −500.5 ± 866j con ζ = 0.5 y ωn = 1000 rad/s; en el osciloscopio el voltaje del capacitor llega a 1.163 V después de cada flanco de subida.
Cambia la resistencia interna de la fuente a 0 o R a 20 Ω y mira cómo se mueven los polos (con R = 20 Ω, ζ ≈ 1: ya no hay sobreimpulso).

Abrir en el simulador →

Ejemplo 3Motor de corriente directa a 12 V

Un motor con

R = 2

Ω,

L = 1

mH,

K = 0.05

V·s/rad,

J = 1 0^{- 4}

kg·m² y

b = 1 0^{- 4}

N·m·s se conecta a una batería de 12 V (0.01 Ω internos). Halla la velocidad final, la corriente final y de arranque, y los polos.

Solución

En régimen permanente las derivadas son cero: $R i + K ω = V$ y $b ω = K i$ . Con $R = 2.01$ Ω: $ω = \frac{K V}{R b + K ^{2}} = \frac{0.05 \cdot 12}{2.01 \cdot 1 0 ^{- 4} + 0.0025} = 222.1$ rad/s $= 2121$ rpm, e $i = bω / K = 0.444$ A.
Al arrancar $ω = 0$ y no hay fuerza contraelectromotriz: la corriente sube hacia $V / R = 5.97$ A, frenada por la inductancia (el pico simulado es 5.83 A). Por eso los motores se protegen contra el arranque.
Polos: $(L s + R) (J s + b) + K^{2} = 0$ → $s^{2} + 2011 s + 27010 = 0$ → $s = - 13.52$ (mecánico, τ = 74 ms) y $s = - 1997$ (eléctrico, 0.5 ms). El lento domina: la velocidad llega al 63 % en unos 74 ms.
Ganancia de voltaje a velocidad: $\frac{2121}{12} = 176.8$ rpm/V.

ω_{\infty} = 2121

rpm,

i_{\infty} = 0.444

A, arranque ≈ 5.8–6 A; polos −13.52 y −1997.

Resuélvelo en Simuveliq

El botón abre un motor M1 con una batería de 12 V. El osciloscopio muestra la velocidad, que llega a unos 2120 rpm a los 0.5 s, y la corriente, con su pico de arranque de 5.83 A que baja a 0.45 A.
Estabilidad: polos −13.52 y −1997 en el subsistema «Motor DC» (el ángulo θ solo acumula y no cuenta). Análisis de control de V1 a la velocidad: G(s) = 4 775 000/(s² + 2011s + 27 010), ganancia en DC 176.8 y valor final del modelo lineal 176.8.

Abrir en el simulador →

Ejemplo 4Calefactor: un sistema térmico de primer orden

Un calefactor de 12 Ω con

C_{t h} = 20

J/K y

R_{t h} = 5

K/W, en un ambiente a 25 °C, se conecta a 12 V. Halla la potencia, la temperatura final, la constante de tiempo y la temperatura a los 100 s.

Solución

Voltaje en el calefactor: $12 \cdot 12/12.01 = 11.990$ V, así que $P = V^{2} / R = 11.98$ W.
Balance de energía: $C_{t h} \dot{T} = P - \frac{T - T _{a}}{R _{t h}}$ . Es un primer orden con $τ = R_{t h} C_{t h} = 100$ s y final $T_{\infty} = T_{a} + P R_{t h} = 25 + 59.9 = 84.9$ °C.
$T (100 s) = 25 + 59.9 (1 - e^{- 1}) = 62.9$ °C.
La potencia depende de $V^{2}$ (no lineal). Linealizando en 12 V: $\frac{d T _{\infty}}{d V} = R_{t h} \frac{2 V _{h}}{R} \cdot \frac{12}{12.01} = 9.98$ °C/V: un volt más subiría unos 10 °C.

P = 11.98

τ = 100

T_{\infty} = 84.9

°C,

T (100 s) = 62.9

°C; ganancia lineal 9.98 °C/V.

Resuélvelo en Simuveliq

El botón abre el calefactor H1 con su batería; la simulación dura 600 s (6τ) y termina en 84.75 °C, a punto de llegar a 84.9. El osciloscopio da la potencia, 11.98 W.
Estabilidad: polo −0.01 con τ = 100 s, «Constante de tiempo térmica τ = C·Rth». Análisis de control de V1 a la temperatura: G(s) = 0.09983/(s + 0.01), ganancia en DC 9.983.

Abrir en el simulador →

Ejemplo 5Tanque: linealizar una salida no lineal

Un tanque de área

A = 0.01

m² descarga por una salida

Q = k h

con

k = 4 \times 1 0^{- 4}

m^2.5/s, y un PID con sensor de presión mantiene el nivel en 0.5 m. ¿Qué constante de tiempo tiene el tanque solo alrededor de ese nivel?

Solución

Balance de volumen: $A \dot{h} = Q_{in} - k h$ . Es no lineal por la raíz.
Cerca de $h_{0}$ se cambia la curva por su tangente: $k h \approx k h_{0} + \frac{k}{2 h _{0}} (h - h_{0})$ . Con $δ h = h - h_{0}$ : $A δ \dot{h} = δ Q_{in} - \frac{k}{2 h _{0}} δ h$ .
Polo: $s = - \frac{k}{2 A h _{0}}$ . En el nivel que alcanza la simulación, $h_{0} = 0.499$ m: $s = - \frac{4 \times 1 0 ^{- 4}}{2 \cdot 0.01 \cdot 0.7064} = - 0.0283$ s⁻¹, $τ = 35.3$ s. A otro nivel el polo cambia: el tanque es más lento cuanto más lleno.
Por eso el modelo lineal y la simulación no coinciden del todo: el escalón de la referencia arranca con el tanque lejos de 0.5 m (sobreimpulso 15.6 % en el modelo lineal y 5.31 % en la simulación).

Alrededor de 0.5 m el tanque tiene un polo en −0.0283 s⁻¹ (

τ = 35.3

s), que depende del nivel.

Resuélvelo en Simuveliq

«Cargar ejemplo» → Control de nivel de un tanque con sensor de presión. En la tabla de subsistemas, «TK1 Tanque — Nivel (k√h linealizado) y válvula» tiene los polos −0.02831 (nivel) y −2 (la válvula, con 0.5 s de respuesta).
En Análisis de control, la tabla compara la simulación no lineal con el modelo lineal: sobreimpulso 5.31 % frente a 15.6 %.

Abrir en el simulador →

Temas relacionados

Sistemas de control: bloques, respuesta y estabilidadLazo abierto y cerrado, función de transferencia, polos y ceros, álgebra de bloques, respuesta de primer y segundo orden, Routh–Hurwitz y error en estado estable.Circuitos eléctricos de corriente directaLey de Ohm, leyes de Kirchhoff, serie y paralelo, transformación Δ–Y, superposición, análisis nodal, Thévenin y transitorios RC y RL.Propiedades termodinámicas y uso de las tablasSustancia pura y fases, campana de saturación, calidad, región de un estado, interpolación simple y doble, líquido comprimido y gas ideal con Pr.

← Todos los temas Ir al simulador