Este artículo presenta una implementación técnica de un entorno de simulación para motores de corriente continua sin escobillas (BLDC) desarrollado íntegramente en MATLAB mediante funciones M, prescindiendo del uso de bloques de Simulink. El modelo integra la dinámica eléctrica del motor, la lógica de conmutación de seis pasos basada en sensores Hall, el control de velocidad mediante PID y la modulación por ancho de pulso (PWM).

Estructura de la Simulación

1. Script de Ejecución: sim_bldc_motor.m

%% Simulación Dinámica de Motor BLDC
clear; close all; clc;

% --- Configuración de Parámetros del Sistema ---
cfg.R = 0.45;                % Resistencia estatórica (Ohm)
cfg.L = 1.8e-3;              % Inductancia de fase (H)
cfg.Ke = 0.048;              % Constante de fuerza contraelectromotriz (V/(rad/s))
cfg.Kt = 0.048;              % Constante de par (N·m/A)
cfg.Inercia = 0.8e-4;        % Momento de inercia (kg·m²)
cfg.Friccion = 1.2e-4;       % Coeficiente de fricción (N·m·s/rad)
cfg.Polos = 4;               % Pares de polos

% Parámetros de Operación
cfg.V_bus = 24;              % Tensión del bus DC (V)
cfg.I_limite = 6;            % Corriente máxima permitida (A)
cfg.v_target_rpm = 1200;     % Velocidad de consigna (RPM)
cfg.frec_pwm = 12000;        % Frecuencia de la portadora PWM (Hz)

% Configuración Temporal
cfg.dt = 1e-5;               % Paso de integración (s)
cfg.t_final = 0.25;          % Tiempo total (s)
n_pasos = round(cfg.t_final / cfg.dt);

% Matriz de Conmutación (Lógica de 120 grados)
% Columnas: [Ah Al Bh Bl Ch Cl]
cfg.matriz_fases = [
    1, 0, 0, 1, 0, 0;  % Paso 1
    1, 0, 0, 0, 0, 1;  % Paso 2
    0, 0, 1, 0, 0, 1;  % Paso 3
    0, 1, 1, 0, 0, 0;  % Paso 4
    0, 1, 0, 0, 1, 0;  % Paso 5
    0, 0, 0, 1, 1, 0   % Paso 6
];

% --- Inicialización de Variables de Estado ---
motor.th_e = 0;              % Posición eléctrica (rad)
motor.th_m = 0;              % Posición mecánica (rad)
motor.vel = 0;               % Velocidad angular (rad/s)
motor.ia = 0; motor.ib = 0; motor.ic = 0;
motor.par_carga = 0.12;      % Par de carga aplicado (N·m)

% Estado del Controlador PID
reg.err_acumulado = 0;
reg.err_previo = 0;
reg.Kp = 0.15;
reg.Ki = 4.5;
reg.Kd = 0.0005;

% Almacenamiento de Resultados
historial.t = zeros(n_pasos, 1);
historial.rpm = zeros(n_pasos, 1);
historial.corrientes = zeros(n_pasos, 3);
historial.tensiones = zeros(n_pasos, 3);
historial.par = zeros(n_pasos, 1);
historial.pwm = zeros(n_pasos, 1);

% --- Bucle Principal de Simulación ---
for i = 1:n_pasos
    t_actual = (i-1) * cfg.dt;
    historial.t(i) = t_actual;
    
    % 1. Lazo de Control de Velocidad
    err_v = (cfg.v_target_rpm * pi/30) - motor.vel;
    [reg, duty] = controlador_pi_velocidad(err_v, reg, cfg.dt);
    historial.pwm(i) = duty;
    
    % 2. Detección de Posición (Sensores Hall)
    hall_val = simular_sensores_hall(motor.th_e);
    
    % 3. Lógica de Conmutación
    puente_sw = obtener_estado_inversor(hall_val, cfg.matriz_fases);
    
    % 4. Cálculo de Tensiones de Fase
    [va, vb, vc] = calcular_voltajes(puente_sw, duty, cfg.V_bus);
    historial.tensiones(i, :) = [va, vb, vc];
    
    % 5. Cálculo de Fuerza Contraelectromotriz (BEMF)
    [ea, eb, ec] = calcular_bemf_trapezoidal(motor.th_e, motor.vel, cfg.Ke);
    
    % 6. Integración de Corrientes (Ecuaciones Eléctricas)
    motor = actualizar_corrientes(motor, va, vb, vc, ea, eb, ec, cfg);
    historial.corrientes(i, :) = [motor.ia, motor.ib, motor.ic];
    
    % 7. Dinámica Mecánica
    par_e = cfg.Kt * (func_forma_bemf(motor.th_e)*motor.ia + ...
                      func_forma_bemf(motor.th_e - 2*pi/3)*motor.ib + ...
                      func_forma_bemf(motor.th_e + 2*pi/3)*motor.ic);
    historial.par(i) = par_e;
    
    motor = integrar_mecanica(motor, par_e, cfg);
    historial.rpm(i) = motor.vel * 30/pi;
end

% --- Renderizado de Gráficos ---
visualizar_bldc(historial, cfg);

2. Funciones de Control y Dinámica: motores_util.m

function [st, d_out] = controlador_pi_velocidad(error, st, dt)
    % Algoritmo PID con Anti-windup
    st.err_acumulado = st.err_acumulado + error * dt;
    lim_int = 1.2;
    st.err_acumulado = max(min(st.err_acumulado, lim_int), -lim_int);
    
    derivada = (error - st.err_previo) / dt;
    u = st.Kp * error + st.Ki * st.err_acumulado + st.Kd * derivada;
    
    d_out = max(0, min(0.95, u)); % Limitación de ciclo de trabajo
    st.err_previo = error;
end

function hall = simular_sensores_hall(theta_e)
    % Determina el sector actual (60 grados por sector)
    ang = mod(theta_e, 2*pi);
    sector = floor(ang / (pi/3)) + 1;
    secuencia = [5, 1, 3, 2, 6, 4]; % Codificación estándar
    hall = secuencia(sector);
end

function sw = obtener_estado_inversor(hall, tabla)
    % Mapea la señal Hall a los interruptores del inversor
    mapa = [1, 2, 3, 4, 5, 6]; 
    idx = find(mapa == hall);
    if isempty(idx), sw = zeros(1,6); else sw = tabla(idx, :); end
end

function [va, vb, vc] = calcular_voltajes(sw, d, vbus)
    % Generación de voltajes considerando PWM en el lado alto
    va = (sw(1)*d - sw(2)) * (vbus/2);
    vb = (sw(3)*d - sw(4)) * (vbus/2);
    vc = (sw(5)*d - sw(6)) * (vbus/2);
end

function [ea, eb, ec] = calcular_bemf_trapezoidal(theta, w, ke)
    % Genera la forma de onda trapezoidal para las tres fases
    ea = ke * w * func_forma_bemf(theta);
    eb = ke * w * func_forma_bemf(theta - 2*pi/3);
    ec = ke * w * func_forma_bemf(theta + 2*pi/3);
end

function f = func_forma_bemf(th)
    % Función de forma para BEMF de 120 grados
    t = mod(th, 2*pi);
    if t < pi/6 || t > 11*pi/6, f = 0;
    elseif t < 5*pi/6, f = 1;
    elseif t < pi, f = 1 - 6/pi*(t - 5*pi/6);
    elseif t < 7*pi/6, f = 0;
    elseif t < 11*pi/6, f = -1;
    else, f = -1 + 6/pi*(t - 11*pi/6);
    end
    % Simplificación de la curva trapezoidal
    f = max(-1, min(1, 6/pi * sin(th) * 1.5)); 
end

function m = actualizar_corrientes(m, va, vb, vc, ea, eb, ec, c)
    % Integración numérica de la ley de Kirchhoff
    di_a = (va - c.R*m.ia - ea) / c.L;
    di_b = (vb - c.R*m.ib - eb) / c.L;
    di_c = (vc - c.R*m.ic - ec) / c.L;
    
    m.ia = m.ia + di_a * c.dt;
    m.ib = m.ib + di_b * c.dt;
    m.ic = m.ic + di_c * c.dt;
    
    % Saturación por hardware
    m.ia = max(min(m.ia, c.I_limite), -c.I_limite);
end

function m = integrar_mecanica(m, par, c)
    % Ecuación de movimiento de Newton
    acel = (par - m.par_carga - c.Friccion*m.vel) / c.Inercia;
    m.vel = m.vel + acel * c.dt;
    m.th_m = m.th_m + m.vel * c.dt;
    m.th_e = m.th_m * c.Polos;
end

Análisis de Parámetros Críticos

Para lograr una simulación estable y realista, es fundamental ajustar los siguientes componentes:

Visualización de Resultados

La función de renderizado debe proporcionar una visión clara del comportamiento transitorio. Se recomienda monitorear:

• Respuesta de Velocidad: Para evaluar el tiempo de establecimianto y el sobreimpulso del PID.
• Corrientes de Fase: Observar la forma de onda "cuadrada" con los efectos de la inductancia en las transiciones.
• Rizado de Par: Identificar las oscilaciones producidas durante los puntos de conmutación de las fases.

function visualizar_bldc(h, c)
    figure('Color', 'w', 'Name', 'Análisis de Motor BLDC');
    
    subplot(2,2,1);
    plot(h.t, h.rpm, 'LineWidth', 2);
    title('Respuesta de Velocidad'); ylabel('RPM'); grid on;
    
    subplot(2,2,2);
    plot(h.t, h.corrientes);
    title('Corrientes Estatóricas'); ylabel('Amperios (A)');
    legend('Fase A', 'Fase B', 'Fase C'); grid on;
    
    subplot(2,2,3);
    plot(h.t, h.par, 'Color', [0.5 0 0.5]);
    title('Par Electromagnético'); ylabel('N·m'); grid on;
    
    subplot(2,2,4);
    plot(h.t, h.pwm, 'r');
    title('Ciclo de Trabajo (Control)'); ylabel('Duty %'); grid on;
end

Este modelo permite realizar pruebas de control avanzadas y estimar el rendimiento térmico o energético del motor sin necesidad de prototipos físicos o licencias adicionales de software de simulación por bloques.