Simulation of PDE nonlinear system using Method of Lines

Hello,

I'm trying to simulate a heat exchanger model that contains PDEs using level 2 s-function. The model gives results but they are not reasonable, I have doubts in the derivative part of the code.

Can anyone help me please, I have been working on this for weeks. I would be very thankful for any kind of help.

function trial_1304(block)
%MSFUNTMPL A Template for a MATLAB S-Function
%   The MATLAB S-function is written as a MATLAB function with the
%   same name as the S-function. Replace 'msfuntmpl' with the name
%   of your S-function.  

%   Copyright 2003-2018 The MathWorks, Inc.
  
%
% The setup method is used to setup the basic attributes of the
% S-function such as ports, parameters, etc. Do not add any other
% calls to the main body of the function.  
%   
setup(block);
  
%endfunction

% Function: setup ===================================================
% Abstract:
%   Set up the S-function block's basic characteristics such as:
%   - Input ports
%   - Output ports
%   - Dialog parameters
%   - Options
% 
%   Required         : Yes
%   C MEX counterpart: mdlInitializeSizes
%
function setup(block)

  % Register the number of ports.
  block.NumInputPorts  = 2;
  block.NumOutputPorts = 3;

  % Register the parameters.
  block.NumDialogPrms     = 11;
  
  % Set up the port properties to be inherited or dynamic.
  block.SetPreCompInpPortInfoToDynamic;
  block.SetPreCompOutPortInfoToDynamic;

  % Override the input port properties.
  block.InputPort(1).Dimensions  = 1;
  block.InputPort(1).DatatypeID  = 0;  % double
  block.InputPort(1).Complexity  = 'Real';
  
  block.InputPort(2).Dimensions  = 1;
  block.InputPort(2).DatatypeID  = 0;  % double
  block.InputPort(2).Complexity  = 'Real';

  % Override the output port properties.  
  block.OutputPort(1).Dimensions  = block.DialogPrm(1).Data;
  block.OutputPort(1).DatatypeID  = 0; % double
  block.OutputPort(1).Complexity  = 'Real';

  block.OutputPort(2).Dimensions  = block.DialogPrm(1).Data;
  block.OutputPort(2).DatatypeID  = 0; % double
  block.OutputPort(2).Complexity  = 'Real';
    
  block.OutputPort(3).Dimensions  = block.DialogPrm(1).Data;
  block.OutputPort(3).DatatypeID  = 0; % double
  block.OutputPort(3).Complexity  = 'Real';

  
  % Set up the continuous states.
  block.NumContStates     = 3*block.DialogPrm(1).Data;

  % Register the sample times.
  %  [0 offset]            : Continuous sample time
  %  [positive_num offset] : Discrete sample time
  %
  %  [-1, 0]               : Inherited sample time
  %  [-2, 0]               : Variable sample time
  block.SampleTimes = [0 0];
  
  % -----------------------------------------------------------------
  % Options
  % -----------------------------------------------------------------
  % Specify if Accelerator should use TLC or call back to the 
  % MATLAB file
  block.SetAccelRunOnTLC(false);
  
  % Specify the block's operating point compliance. The block operating 
  % point is used during the containing model's operating point save/restore)
  % The allowed values are:
  %   'Default' : Same the block's operating point as of a built-in block
  %   'UseEmpty': No data to save/restore in the block operating point
  %   'Custom'  : Has custom methods for operating point save/restore
  %                 (see GetOperatingPoint/SetOperatingPoint below)
  %   'Disallow': Error out when saving or restoring the block operating point.
  block.OperatingPointCompliance = 'Default';
  
  % -----------------------------------------------------------------
  % The MATLAB S-function uses an internal registry for all
  % block methods. You should register all relevant methods
  % (optional and required) as illustrated below. You may choose
  % any suitable name for the methods and implement these methods
  % as local functions within the same file.
  % -----------------------------------------------------------------

  %
  % SetInputPortSamplingMode:
  %   Functionality    : Check and set input and output port 
  %                      attributes and specify whether the port is operating 
  %                      in sample-based or frame-based mode
  %   C MEX counterpart: mdlSetInputPortFrameData.
  %   (The DSP System Toolbox is required to set a port as frame-based)
  %
block.RegBlockMethod('SetInputPortSamplingMode',@SetInputPortSamplingMode);
 
  %
  % PostPropagationSetup:
  %   Functionality    : Set up the work areas and the state variables. You can
  %                      also register run-time methods here.
  %   C MEX counterpart: mdlSetWorkWidths
  %
  block.RegBlockMethod('PostPropagationSetup', @DoPostPropSetup);

  % -----------------------------------------------------------------
  % Register methods called at run-time
  % -----------------------------------------------------------------
  
  % 
  % ProcessParameters:
  %   Functionality    : Call to allow an update of run-time parameters.
  %   C MEX counterpart: mdlProcessParameters
  %  
  block.RegBlockMethod('ProcessParameters', @ProcessPrms);

  % 
  % InitializeConditions:
  %   Functionality    : Call to initialize the state and the work
  %                      area values.
  %   C MEX counterpart: mdlInitializeConditions
  % 
  block.RegBlockMethod('InitializeConditions', @InitializeConditions);
  
  % 
  % Start:
  %   Functionality    : Call to initialize the state and the work
  %                      area values.
  %   C MEX counterpart: mdlStart
  %
  block.RegBlockMethod('Start', @Start);

  % 
  % Outputs:
  %   Functionality    : Call to generate the block outputs during a
  %                      simulation step.
  %   C MEX counterpart: mdlOutputs
  %
  block.RegBlockMethod('Outputs', @Outputs);

  % 
  % Update:
  %   Functionality    : Call to update the discrete states
  %                      during a simulation step.
  %   C MEX counterpart: mdlUpdate
  %
  block.RegBlockMethod('Update', @Update);

  % 
  % Derivatives:
  %   Functionality    : Call to update the derivatives of the
  %                      continuous states during a simulation step.
  %   C MEX counterpart: mdlDerivatives
  %
  block.RegBlockMethod('Derivatives', @Derivatives);
  
 
  % 
  % SimStatusChange:
  %   Functionality    : Call when simulation enters pause mode
  %                      or leaves pause mode.
  %   C MEX counterpart: mdlSimStatusChange
  %
  block.RegBlockMethod('SimStatusChange', @SimStatusChange);
  
  % 
  % Terminate:
  %   Functionality    : Call at the end of a simulation for cleanup.
  %   C MEX counterpart: mdlTerminate
  %
  block.RegBlockMethod('Terminate', @Terminate);

  %
  % GetOperatingPoint:
  %   Functionality    : Return the operating point of the block.
  %   C MEX counterpart: mdlGetOperatingPoint
  %
  block.RegBlockMethod('GetOperatingPoint', @GetOperatingPoint);
  
  %
  % SetOperatingPoint:
  %   Functionality    : Set the operating point data of the block using
  %                       from the given value.
  %   C MEX counterpart: mdlSetOperatingPoint
  %
  block.RegBlockMethod('SetOperatingPoint', @SetOperatingPoint);

  % -----------------------------------------------------------------
  % Register the methods called during code generation.
  % -----------------------------------------------------------------
  
  %
  % WriteRTW:
  %   Functionality    : Write specific information to model.rtw file.
  %   C MEX counterpart: mdlRTW
  %
  block.RegBlockMethod('WriteRTW', @WriteRTW);
%endfunction

% -------------------------------------------------------------------
% The local functions below are provided to illustrate how you may implement
% the various block methods listed above.
% -------------------------------------------------------------------



function ProcessPrms(block)

  block.AutoUpdateRuntimePrms;
 
%endfunction
function SetInputPortSamplingMode(block, port, sp)
    block.InputPort(port).SamplingMode  = sp; 
    block.OutputPort(1).SamplingMode    = sp; 
    block.OutputPort(2).SamplingMode    = sp;
    block.OutputPort(3).SamplingMode    = sp;

        
function DoPostPropSetup(block)
  block.NumDworks = 3;
  
  block.Dwork(1).Name            = 'Tt';
  block.Dwork(1).Dimensions      = block.DialogPrm(1).Data;
  block.Dwork(1).DatatypeID      = 0;      % double
  block.Dwork(1).Complexity      = 'Real'; % real
  block.Dwork(1).UsedAsDiscState = true;
  
  block.Dwork(2).Name            = 'Ts';
  block.Dwork(2).Dimensions      = block.DialogPrm(1).Data;
  block.Dwork(2).DatatypeID      = 0;      % double
  block.Dwork(2).Complexity      = 'Real'; % real
  block.Dwork(2).UsedAsDiscState = true;
  
  block.Dwork(3).Name            = 'Tm';
  block.Dwork(3).Dimensions      = block.DialogPrm(1).Data;
  block.Dwork(3).DatatypeID      = 0;      % double
  block.Dwork(3).Complexity      = 'Real'; % real
  block.Dwork(3).UsedAsDiscState = true;
  % Register all tunable parameters as runtime parameters.
  block.AutoRegRuntimePrms;

%endfunction

function InitializeConditions(block)

%endfunction

function Start(block)

n = block.DialogPrm(1).Data;

stateSize = 3*n;
block.ContStates.Data = zeros(stateSize,1);

block.Dwork(1).Data(1) = block.DialogPrm(3).Data ; 
block.Dwork(2).Data(1) = block.DialogPrm(6).Data ; 
block.Dwork(3).Data(1) = block.DialogPrm(9).Data ; 
for i = 2:n-1
    block.Dwork(1).Data(i) = block.DialogPrm(4).Data ; 
    block.Dwork(2).Data(i) = block.DialogPrm(7).Data ; 
    block.Dwork(3).Data(i) = block.DialogPrm(10).Data ; 
end
block.Dwork(1).Data(n) = block.DialogPrm(5).Data ; 
block.Dwork(2).Data(n) = block.DialogPrm(8).Data ; 
block.Dwork(3).Data(n) = block.DialogPrm(11).Data ; 


block.ContStates.Data(1:n)          = block.Dwork(1).Data;
block.ContStates.Data(n+1:2*n)      = block.Dwork(2).Data;
block.ContStates.Data(2*n+1:3*n)    = block.Dwork(3).Data;

%endfunction


function Outputs(block)
n = block.DialogPrm(1).Data;

  
  block.OutputPort(1).Data = block.ContStates.Data(1:n);
  block.OutputPort(2).Data = block.ContStates.Data(n+1:2*n);
  block.OutputPort(3).Data = block.ContStates.Data(2*n+1:3*n);

  
%endfunction

function Update(block)
  
    block.Dwork(1).Data(1) = block.InputPort(1).Data;
    block.Dwork(2).Data(1) = block.InputPort(2).Data;
  
%endfunction

function Derivatives(block)

    c = 2;
    Ut = 0.284;
    Us = 0.284;
    AV = 157;
    rhom = 7700;
    cm = 0.5;
    thm = 3.048e-3;
    vt = 0.5;
    vs = 0.5;
    rhot = 0.5;
    rhos = 0.5;
    

    L = block.DialogPrm(2).Data;
    n = block.DialogPrm(1).Data;
    z = linspace(0,L,n); 

%     Tt = block.ContStates.Data(1:n);
%     Ts = block.ContStates.Data(n+1:2*n);
%     Tm = block.ContStates.Data(2*n+1:3*n);  

%     dTtdt = zeros(n,1);
%     dTsdt = zeros(n,1);
%     dTmdt = zeros(n,1);

    dTtdz = zeros(n,1);
    dTsdz = zeros(n,1);
    dTmdz = zeros(n,1);
   
    for i=2:n
        dTtdz(i) = ( block.Dwork(1).Data(i)- block.Dwork(1).Data(i-1))/(z(i)-z(i-1));
        dTsdz(i) = ( block.Dwork(2).Data(i)- block.Dwork(2).Data(i-1))/(z(i)-z(i-1));
        dTmdz(i) = ( block.Dwork(3).Data(i)- block.Dwork(3).Data(i-1))/(z(i)-z(i-1));
    end 
%         
%     for i=2:n
%         dTtdt(i) = -vt*dTtdz(i)+Ut*AV*(block.Dwork(3).Data(i)-block.Dwork(1).Data(i))/(rhot*c);
%         dTsdt(i) = vs*dTsdz(i)-Us*AV*(block.Dwork(2).Data(i)-block.Dwork(3).Data(i))/(rhos*c);
%         dTmdt(i) = -(Us*(block.Dwork(2).Data(i)-block.Dwork(3).Data(i))-Ut*(block.Dwork(3).Data(i)-block.Dwork(1).Data(i)))/(rhom*cm*thm);
%     end


      block.Derivatives.Data(1)     = 0;
      block.Derivatives.Data(n+1)   = 0;
      block.Derivatives.Data(2*n+1) = 0;
 
  for i = 2:n
      block.Derivatives.Data(i)     = -vt*dTtdz(i)+Ut*AV*(block.Dwork(3).Data(i)-block.Dwork(1).Data(i))/(rhot*c);
      block.Derivatives.Data(n+i)   = vs*dTsdz(i)-Us*AV*(block.Dwork(2).Data(i)-block.Dwork(3).Data(i))/(rhos*c);
      block.Derivatives.Data(2*n+i) = -(Us*(block.Dwork(2).Data(i)-block.Dwork(3).Data(i))-Ut*(block.Dwork(3).Data(i)-block.Dwork(1).Data(i)))/(rhom*cm*thm);
  end 

%endfunction

    
function Terminate(block)

disp(['Terminating the block with handle ' num2str(block.BlockHandle) '.']);

%endfunction