Code covered by the BSD License  

Highlights from
Discrete Lorenz Water Wheel Simulation

image thumbnail

Discrete Lorenz Water Wheel Simulation

by

 

This is a simulation of an 8 tank lorenze water wheel. The water speed is manually controlled.

lorenze_dynamics(t,x,u,flag,P)
function [sys,x0,str,ts,simStateCompliance] = lorenze_dynamics(t,x,u,flag,P)



switch flag,

  %%%%%%%%%%%%%%%%%%
  % Initialization %
  %%%%%%%%%%%%%%%%%%
  case 0,
    [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes(P);

  %%%%%%%%%%%%%%%
  % Derivatives %
  %%%%%%%%%%%%%%%
  case 1,
    sys=mdlDerivatives(t,x,u,P);

  %%%%%%%%%%
  % Update %
  %%%%%%%%%%
  case 2,
    sys=mdlUpdate(t,x,u);

  %%%%%%%%%%%
  % Outputs %
  %%%%%%%%%%%
  case 3,
    sys=mdlOutputs(t,x,u,P);

  %%%%%%%%%%%%%%%%%%%%%%%
  % GetTimeOfNextVarHit %
  %%%%%%%%%%%%%%%%%%%%%%%
  case 4,
    sys=mdlGetTimeOfNextVarHit(t,x,u);

  %%%%%%%%%%%%%
  % Terminate %
  %%%%%%%%%%%%%
  case 9,
    sys=mdlTerminate(t,x,u);

  %%%%%%%%%%%%%%%%%%%%
  % Unexpected flags %
  %%%%%%%%%%%%%%%%%%%%
  otherwise
    DAStudio.error('Simulink:blocks:unhandledFlag', num2str(flag));

end

% end sfuntmpl

%
%=============================================================================
% mdlInitializeSizes
% Return the sizes, initial conditions, and sample times for the S-function.
%=============================================================================
%
function [sys,x0,str,ts,simStateCompliance]=mdlInitializeSizes(P)


%
% call simsizes for a sizes structure, fill it in and convert it to a
% sizes array.
%
% Note that in this example, the values are hard coded.  This is not a
% recommended practice as the characteristics of the block are typically
% defined by the S-function parameters.
%
sizes = simsizes;

sizes.NumContStates  = 17;
sizes.NumDiscStates  = 0;
sizes.NumOutputs     = 17;
sizes.NumInputs      = 18;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;   % at least one sample time is needed

sys = simsizes(sizes);

%
% initialize the initial conditions
%
x0_theta  = linspace(0,(2*pi - 2*pi/8),8)';
x0_heights = P.TH*.5*ones(8,1);
x0_velWheel = 0;
x0 = [x0_theta;x0_heights;x0_velWheel];

%
% str is always an empty matrix
%
str = [];

%
% initialize the array of sample times
%
ts  = [0 0];

% Specify the block simStateCompliance. The allowed values are:
%    'UnknownSimState', < The default setting; warn and assume DefaultSimState
%    'DefaultSimState', < Same sim state as a built-in block
%    'HasNoSimState',   < No sim state
%    'DisallowSimState' < Error out when saving or restoring the model sim state
simStateCompliance = 'UnknownSimState';

% end mdlInitializeSizes

%==========================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u,P)


 % positions 
 
 theta_vector(:,1) = u(1:8,1);
 
 r_vector = ones(8,1)*P.rmax;
 % water heights
 
 h_vector(:,1) = u(9:16,1);
 

 
 % Velocity of the wheel
 
 v_wheel = u(17);
 
 
 density = 1000;     % kg/m^3
 Ac = pi*(P.TW/2)^2;
 mass = Ac.*h_vector*density + .00005;
 
 alpha = AlphaCalc(r_vector,theta_vector,mass,P);  % Angular Acceleration 
 
 v_dot = alpha - P.b*v_wheel;
 

  % Velocity of the water 
 
 gamma = u(18);    % Valve position
 
 VolFlow = P.VolFlowMax*gamma;
 
 angle_check = abs(asin(P.TW/2/P.rmax));
  
 upperBound = 2*pi - angle_check;
 lowerBound = angle_check;
 h_dot = h_vector.^0;
 
 for i = 1:length(h_vector)
     
     if theta_vector(i)>= upperBound || theta_vector(i) <= lowerBound
         
         h_dot(i) = VolFlow/Ac;
         
     else
         h_dot(i) = 0;
         
     end
 end
 
h_out = (P.outFlowConst)/Ac.*h_vector.^2;

h_dot = h_dot - h_out;
theta_dot = ones(8,1)*v_wheel;
  
 
      
  
sys = [theta_dot; h_dot;  v_dot];

% end mdlDerivatives

%
%=============================================================================
% mdlUpdate
% Handle discrete state updates, sample time hits, and major time step
% requirements.
%=============================================================================
%
function sys=mdlUpdate(t,x,u)

sys = [];

% end mdlUpdate

%
%=============================================================================
% mdlOutputs
% Return the block outputs.
%=============================================================================
%
function sys=mdlOutputs(t,x,u,P)

  theta_vector(:,1) = mod(x(1:8,1),2*pi);
   h_vector(:,1) = x(9:16,1);
   v_wheel = x(17);
      
  
sys = [theta_vector; h_vector; v_wheel];

% end mdlOutputs

%
%=============================================================================
% mdlGetTimeOfNextVarHit
% Return the time of the next hit for this block.  Note that the result is
% absolute time.  Note that this function is only used when you specify a
% variable discrete-time sample time [-2 0] in the sample time array in
% mdlInitializeSizes.
%=============================================================================
%
function sys=mdlGetTimeOfNextVarHit(t,x,u)

sampleTime = 1;    %  Example, set the next hit to be one second later.
sys = t + sampleTime;

% end mdlGetTimeOfNextVarHit

%
%=============================================================================
% mdlTerminate
% Perform any end of simulation tasks.
%=============================================================================
%
function sys=mdlTerminate(t,x,u)

sys = [];

% end mdlTerminate

Contact us