представляет тепловую модель в Matlab с использованием ODE - PullRequest
0 голосов
/ 06 октября 2018

Я работаю над симуляцией переменного тока так же, как показано в этой ссылке .Мне удалось преобразовать его в переменный ток, перевернув некоторые уравнения и повторно инициализировав некоторые значения.Теперь мне нужно сделать это в Matlab, чтобы моя система работала плавно и быстро, без необходимости запускать simulink для каждой итерации.Я прочитал ODE45 и решил использовать его, во-первых, я преобразовал модель на веб-странице в функцию следующим образом:

    function dTroom_dt = odes_thermal(Troom,Tout)
    TinIC = 26;
    r2d = 180/pi;
    % -------------------------------
    % Define the house geometry
    % -------------------------------
    % House length = 30 m
    lenHouse = 30;
    % House width = 10 m
    widHouse = 10;
    % House height = 4 m
    htHouse = 4;
    % Roof pitch = 40 deg
    pitRoof = 40/r2d;
    % Number of windows = 6
    numWindows = 6;
    % Height of windows = 1 m
    htWindows = 1;
    % Width of windows = 1 m
    widWindows = 1;
    windowArea = numWindows*htWindows*widWindows;
    wallArea = 2*lenHouse*htHouse + 2*widHouse*htHouse + ...
        2*(1/cos(pitRoof/2))*widHouse*lenHouse + ...
        tan(pitRoof)*widHouse - windowArea;
    % -------------------------------
    % Define the type of insulation used
    % -------------------------------
    % Glass wool in the walls, 0.2 m thick
    % k is in units of J/sec/m/C - convert to J/hr/m/C multiplying by 3600
    kWall = 0.038*3600;   % hour is the time unit
    LWall = .2;
    RWall = LWall/(kWall*wallArea);
    % Glass windows, 0.01 m thick
    kWindow = 0.78*3600;  % hour is the time unit
    LWindow = .01;
    RWindow = LWindow/(kWindow*windowArea);
    % -------------------------------
    % Determine the equivalent thermal resistance for the whole building
    % -------------------------------
    Req = RWall*RWindow/(RWall + RWindow);
    % c = cp of air (273 K) = 1005.4 J/kg-K
    c = 1005.4;
    % -------------------------------
    % Enter the temperature of the heated air
    % -------------------------------
    % The air exiting the heater has a constant temperature which is a heater
    % property. THeater =20 deg C
    THeater = 20;
    % Air flow rate Mdot = 1 kg/sec = 3600 kg/hr
    Mdot = 3600;  % hour is the time unit
    % -------------------------------
    % Determine total internal air mass = M
    % -------------------------------
    % Density of air at sea level = 1.2250 kg/m^3
    densAir = 1.2250;
    M = (lenHouse*widHouse*htHouse+tan(pitRoof)*widHouse*lenHouse)*densAir;
    %%

    dQlosses_dt= (THeater-Troom)*Mdot*c;

    dQlosses_dt=(Troom-Tout)/Req;

    dTroom_dt=(1/(M*c))*(dQlosses_dt);

    end

Затем я просто вызвал ее, используя следующие строки:

    clear;
    t=1:1:1440;
    Troom_0=25;
    [t,Troom] = ode45('odes_thermal',[t(1) t(3)],Troom_0);
    figure(1);
    plot(t,Troom);

когда я запускаю его, вот что я получаю: enter image description here

Мой вопрос заключается в том, как включить в него состояние или отключить функции.В Simulink они использовали реле, которое будет выдавать сигнал состояния, который будет 1 или 0 в зависимости от температуры, вот как они это сделали:

enter image description here

Я знаю, что это может быть представлено как "dQlosses_dt = (THeater-Troom) Mdot c * status;"но я не знаю, как обновить состояние для каждого временного образца, так как нет возможности проверить Troom в каждом временном интервале. Можете ли вы мне помочь в этом или хотя бы с чего начать?

Редактировать:Следуя комментарию Pacta, я попытался превратить его в многовариантную систему с двумя выходами Troom и Qlosses.Вот новые функции, но я получаю сообщение об ошибке (слишком много входных аргументов в odes_thermal2)

    clear;
    t=1:1:1440;
    Troom_0=25;Qlosses_0=10;
    X0=[Troom_0;Qlosses_0];
    [t,X] = ode45('odes_thermal2',[t(1) t(end)],X0);
    Troom=X(1);
    Qlosses=X(2);
    figure(1);
    plot(t,Troom);
    ylabel('Troom');
    xlabel('t');

И функция ODE имеет следующие изменения:

    function dX_dt = odes_thermal2(X)
    TinIC = 26;
    r2d = 180/pi;

и

    X(2)= (THeater-Troom)*Mdot*c;

    dQlosses_dt=(Troom-Tout)/Req;

    X(1)=(1/(M*c))*(dQlosses_dt);

Спасибо

...