Simulation of PDE nonlinear system using Method of Lines

1 次查看（过去 30 天）
Mohammed Hamza
Mohammed Hamza 2022-4-12
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

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 PDE Solvers 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by